MRS Longitudinal Modelling
Meaghan Perdue
12/14/2021
Organize data
Merge output data from multiple batches, select metabolite columns, add participant info and other data
acg_batch1 <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/ACG_LCM_CorrectedTable_12-14-2021-16-15.csv")
acg_batch2 <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/ACG_LCM_CorrectedTable_12-14-2021-16-21.csv")
acg_LCM_corr <- bind_rows(acg_batch1, acg_batch2)
anyDuplicated(acg_LCM_corr$File_ID) #check for duplicated values## [1] 106acg_LCM_corr<-unique(acg_LCM_corr) #remove duplicate rows
acg_LCM_corr["4", "Cho_GPC_PCh_conc_molar"]<-NA #Cho conc invalid for PS10008_PS14_046_2386 due to spike in spectrum
acg_LCM_corr["4", "Cho_GPC_PCh_conc_molal"]<-NA #Cho conc invalid for PS10008_PS14_046_2386 due to spike in spectrum
write.csv(acg_LCM_corr, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/ACG_LCM_CorrectedTable_12-14-2021_all.csv", row.names = FALSE)
acg_molal <- acg_LCM_corr %>%
dplyr::select(File_ID, fGM, fWM, fCSF, NAA_conc_molal, Cre_PCr_conc_molal, Cho_GPC_PCh_conc_molal, Ins_conc_molal, Glu_conc_molal, Glu_Gln_conc_molal)
acg_demo <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/acg_analysis_subs.csv")
acg_dat <- left_join(acg_demo, acg_molal, by = "File_ID")
write.csv(acg_dat, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/acg_data_molal_n113.csv", row.names = FALSE)
lag_batch1 <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/LAG_LCM_CorrectedTable_12-13-2021-15-18.csv")
lag_batch2 <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/LAG_LCM_CorrectedTable_12-14-2021-11-41.csv")
lag_batch3 <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/LAG_LCM_CorrectedTable_12-14-2021-12-18.csv")
lag_LCM_corr <- bind_rows(lag_batch1, lag_batch2, lag_batch3)
write.csv(lag_LCM_corr, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/LAG_LCM_CorrectedTable_12-14-2021_all.csv", row.names = FALSE)
lag_molal <- lag_LCM_corr %>%
dplyr::select(File_ID, fGM, fWM, fCSF, NAA_conc_molal, Cre_PCr_conc_molal, Cho_GPC_PCh_conc_molal, Ins_conc_molal, Glu_conc_molal, Glu_Gln_conc_molal)
lag_demo <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/lag_analysis_subs.csv")
lag_dat <- left_join(lag_demo, lag_molal, by = "File_ID")
lag_dat$subj_id<-as.character(lag_dat$subj_id) #make subj_id character type
write.csv(lag_dat, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/lag_data_molal_n321.csv", row.names = FALSE)##Check data distribution & outliers Create plots to check data distribution and identify outliers
ACG
plot_age <- ggplot(acg_dat, aes(mri_age_y))
plot_age + geom_density()
#age skewed
plot_naa <- ggplot(acg_dat, aes(NAA_conc_molal))
plot_naa + geom_density()
plot_naa + geom_boxplot()
rosnerTest(acg_dat$NAA_conc_molal, k = 4)## $distribution
## [1] "Normal"
##
## $statistic
## R.1 R.2 R.3 R.4
## 2.759147 2.825357 2.694361 2.721635
##
## $sample.size
## [1] 113
##
## $parameters
## k
## 4
##
## $alpha
## [1] 0.05
##
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4
## 3.425263 3.422302 3.419309 3.416284
##
## $n.outliers
## [1] 0
##
## $alternative
## [1] "Up to 4 observations are not\n from the same Distribution."
##
## $method
## [1] "Rosner's Test for Outliers"
##
## $data
## [1] 13.36240 15.05227 14.57780 13.83339 14.71206 16.20111 15.44740 16.98446
## [9] 15.67072 15.95324 14.87793 12.74216 16.29728 14.71133 14.66048 13.28760
## [17] 14.59532 13.55773 16.85293 15.33699 12.45017 15.04043 15.27049 11.24177
## [25] 14.98401 14.18774 15.83415 13.91788 17.10858 12.09629 14.70242 14.40906
## [33] 16.46456 15.58993 15.03289 14.84971 11.30097 13.87568 15.24533 15.99486
## [41] 14.53574 12.16795 14.20751 14.18561 16.47436 14.52629 15.00648 16.20561
## [49] 13.22605 13.57437 14.05939 13.81423 17.33147 16.31303 16.29984 16.85100
## [57] 15.24641 14.93994 15.14031 15.13625 18.29604 13.89531 13.42919 15.16179
## [65] 13.87656 15.48760 12.64660 16.01636 13.64857 16.43270 15.59565 11.64927
## [73] 15.46789 14.18371 12.72496 16.73666 16.33378 15.15593 15.86369 15.28939
## [81] 15.31529 14.67805 13.38081 15.01442 13.01089 17.11696 14.54570 15.04389
## [89] 14.64002 15.57552 15.21758 15.80866 15.16132 14.92979 14.18939 16.96225
## [97] 13.86479 15.05179 15.22074 15.65731 15.51013 14.39392 15.03669 16.39137
## [105] 14.47850 14.65614 15.94476 16.78988 16.07591 13.54665 13.48426 14.73463
## [113] 15.68117
##
## $data.name
## [1] "acg_dat$NAA_conc_molal"
##
## $bad.obs
## [1] 0
##
## $all.stats
## i Mean.i SD.i Value Obs.Num R.i+1 lambda.i+1 Outlier
## 1 0 14.88986 1.322180 11.24177 24 2.759147 3.425263 FALSE
## 2 1 14.92243 1.281772 11.30097 37 2.825357 3.422302 FALSE
## 3 2 14.95506 1.239990 18.29604 61 2.694361 3.419309 FALSE
## 4 3 14.92469 1.203472 11.64927 72 2.721635 3.416284 FALSE
##
## attr(,"class")
## [1] "gofOutlier"plot_cre <- ggplot(acg_dat, aes(Cre_PCr_conc_molal))
plot_cre + geom_density()
plot_cre + geom_boxplot()
rosnerTest(acg_dat$Cre_PCr_conc_molal, k = 5)## $distribution
## [1] "Normal"
##
## $statistic
## R.1 R.2 R.3 R.4 R.5
## 3.394290 2.833607 2.707716 2.578893 2.649738
##
## $sample.size
## [1] 113
##
## $parameters
## k
## 5
##
## $alpha
## [1] 0.05
##
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 lambda.5
## 3.425263 3.422302 3.419309 3.416284 3.413225
##
## $n.outliers
## [1] 0
##
## $alternative
## [1] "Up to 5 observations are not\n from the same Distribution."
##
## $method
## [1] "Rosner's Test for Outliers"
##
## $data
## [1] 10.869437 10.445542 14.138594 11.979460 11.842544 12.567296 12.183610
## [8] 13.315617 13.532807 13.910100 12.102060 9.609688 12.566391 11.689508
## [15] 12.650725 10.976876 11.888348 11.823341 12.975541 12.087174 10.745566
## [22] 12.005758 13.958950 11.012768 13.437536 12.359373 12.418194 10.775703
## [29] 13.584172 11.251205 11.993304 12.190118 14.756858 13.780875 12.273252
## [36] 11.243397 8.406874 11.279114 11.602557 13.761404 12.044563 10.276615
## [43] 10.918920 10.896575 11.867174 12.189638 12.580644 12.364280 11.517322
## [50] 10.442494 11.903800 12.566843 12.411654 12.733123 12.852863 14.066514
## [57] 10.742053 11.977272 12.534036 12.467809 13.778825 11.508522 11.740032
## [64] 11.569008 10.891687 12.474281 11.395342 12.716648 12.033736 13.646045
## [71] 12.873670 9.632426 11.901341 12.610355 9.187364 13.087309 12.160354
## [78] 11.717870 12.392428 11.411606 12.417766 11.214681 11.894967 11.689931
## [85] 11.069854 13.648131 12.837410 11.764882 12.329981 11.666319 12.825568
## [92] 13.032058 12.630243 12.355787 11.038718 9.731136 12.044142 11.736949
## [99] 11.772225 11.427900 12.588592 11.074387 11.520514 12.986814 12.563423
## [106] 12.334515 13.095183 12.401466 11.143206 11.414014 11.499070 12.331493
## [113] 11.900989
##
## $data.name
## [1] "acg_dat$Cre_PCr_conc_molal"
##
## $bad.obs
## [1] 0
##
## $all.stats
## i Mean.i SD.i Value Obs.Num R.i+1 lambda.i+1 Outlier
## 1 0 12.03590 1.0691568 8.406874 37 3.394290 3.425263 FALSE
## 2 1 12.06830 1.0167044 9.187364 75 2.833607 3.422302 FALSE
## 3 2 12.09426 0.9833373 14.756858 33 2.707716 3.419309 FALSE
## 4 3 12.07005 0.9540396 9.609688 12 2.578893 3.416284 FALSE
## 5 4 12.09263 0.9284692 9.632426 72 2.649738 3.413225 FALSE
##
## attr(,"class")
## [1] "gofOutlier"plot_cho <- ggplot(acg_dat, aes(Cho_GPC_PCh_conc_molal))
plot_cho + geom_density()## Warning: Removed 1 rows containing non-finite values (stat_density).
plot_cho + geom_boxplot()## Warning: Removed 1 rows containing non-finite values (stat_boxplot).
rosnerTest(acg_dat$Cho_GPC_PCh_conc_molal, k = 4)## Warning in rosnerTest(acg_dat$Cho_GPC_PCh_conc_molal, k = 4): 1 observations
## with NA/NaN/Inf in 'x' removed.## $distribution
## [1] "Normal"
##
## $statistic
## R.1 R.2 R.3 R.4
## 2.692059 2.486294 2.542571 2.587426
##
## $sample.size
## [1] 112
##
## $parameters
## k
## 4
##
## $alpha
## [1] 0.05
##
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4
## 3.422302 3.419309 3.416284 3.413225
##
## $n.outliers
## [1] 0
##
## $alternative
## [1] "Up to 4 observations are not\n from the same Distribution."
##
## $method
## [1] "Rosner's Test for Outliers"
##
## $data
## [1] 3.054728 2.338141 2.483788 2.615586 2.904744 3.029351 2.945605 2.946224
## [9] 2.469310 2.647417 2.329234 2.525650 2.964839 2.377024 2.316960 2.388015
## [17] 2.587590 2.964150 2.987201 1.884924 2.294818 2.160623 1.723047 2.273437
## [25] 2.486538 2.364455 2.110030 3.042283 2.145223 2.626855 2.657139 3.307897
## [33] 2.993478 2.528009 2.694001 2.751861 2.169544 2.576529 2.777803 2.298135
## [41] 2.093485 2.245856 1.887058 3.332061 2.683945 2.960541 2.306091 2.400287
## [49] 2.262987 2.373107 2.418340 2.428692 2.651362 2.724707 2.627568 3.213338
## [57] 2.541659 2.436647 2.503022 3.168533 2.273151 2.230455 2.575725 2.463493
## [65] 2.964749 2.427496 2.382196 2.398624 2.592440 2.075228 1.967159 2.210616
## [73] 2.178031 1.572637 2.265866 2.193475 2.435427 2.749229 2.738993 2.096637
## [81] 2.376787 2.150157 2.031053 1.736687 2.549619 2.032952 2.153716 2.354125
## [89] 2.585598 2.624180 2.597158 2.503254 2.679817 2.109408 2.755517 2.618999
## [97] 2.254582 2.379855 2.105926 2.389925 2.621861 2.349337 2.962557 2.620008
## [105] 2.643999 2.559238 3.322683 2.403398 2.301871 2.505683 2.762249 2.704160
##
## $data.name
## [1] "acg_dat$Cho_GPC_PCh_conc_molal"
##
## $bad.obs
## [1] 1
##
## $all.stats
## i Mean.i SD.i Value Obs.Num R.i+1 lambda.i+1 Outlier
## 1 0 2.495888 0.3429536 1.572637 75 2.692059 3.422302 FALSE
## 2 1 2.504206 0.3329678 3.332061 45 2.486294 3.419309 FALSE
## 3 2 2.496680 0.3248691 3.322683 108 2.542571 3.416284 FALSE
## 4 3 2.489102 0.3164519 3.307897 33 2.587426 3.413225 FALSE
##
## attr(,"class")
## [1] "gofOutlier"plot_ins <- ggplot(acg_dat, aes(Ins_conc_molal))
plot_ins + geom_density()
shapiro.test(acg_dat$Ins_conc_molal) #not normal##
## Shapiro-Wilk normality test
##
## data: acg_dat$Ins_conc_molal
## W = 0.96211, p-value = 0.002754plot_ins + geom_boxplot()
rosnerTest(acg_dat$Ins_conc_molal, k = 3) #outlier value 1.471244, observation 84 (PS18_013) - spectrum showed issue, exclude from analysis ## $distribution
## [1] "Normal"
##
## $statistic
## R.1 R.2 R.3
## 3.753416 3.044570 3.051178
##
## $sample.size
## [1] 113
##
## $parameters
## k
## 3
##
## $alpha
## [1] 0.05
##
## $crit.value
## lambda.1 lambda.2 lambda.3
## 3.425263 3.422302 3.419309
##
## $n.outliers
## [1] 1
##
## $alternative
## [1] "Up to 3 observations are not\n from the same Distribution."
##
## $method
## [1] "Rosner's Test for Outliers"
##
## $data
## [1] 7.633881 8.178435 7.218240 7.696923 8.696692 8.600324 9.402207
## [8] 8.875539 6.779586 5.935284 10.304907 5.028532 7.058793 9.067489
## [15] 6.700483 7.361395 8.124760 8.682546 6.553380 6.320046 7.315123
## [22] 7.446418 7.281237 7.102747 9.505161 7.839894 7.887706 7.616099
## [29] 7.601724 6.347168 7.912922 5.656884 9.656654 9.672370 8.341248
## [36] 8.816154 5.196891 6.304340 6.789199 8.018694 7.594973 7.416233
## [43] 7.406756 3.198494 2.985922 9.898020 8.079535 8.703354 7.157223
## [50] 6.213808 8.399767 9.174211 9.656691 9.335033 10.080614 4.813291
## [57] 9.728369 7.521052 7.441993 8.051283 10.887902 8.695335 5.872553
## [64] 6.418561 5.242212 8.749670 7.765345 8.185832 7.950928 8.394652
## [71] 9.583041 6.202420 7.829477 6.209769 4.285017 8.677590 7.775045
## [78] 9.334317 8.684385 7.236647 7.742767 6.219638 9.105159 1.471244
## [85] 5.218826 9.407651 7.634559 6.064850 8.785683 3.770173 7.618745
## [92] 8.763466 7.123140 6.399259 6.804541 5.855502 9.364210 6.203701
## [99] 10.565682 9.856594 7.453487 7.728481 7.846516 9.947994 8.640216
## [106] 7.758280 8.476271 8.672477 6.437713 5.739577 6.673856 7.596796
## [113] 7.229168
##
## $data.name
## [1] "acg_dat$Ins_conc_molal"
##
## $bad.obs
## [1] 0
##
## $all.stats
## i Mean.i SD.i Value Obs.Num R.i+1 lambda.i+1 Outlier
## 1 0 7.588864 1.629881 1.471244 84 3.753416 3.425263 TRUE
## 2 1 7.643485 1.529793 2.985922 45 3.044570 3.422302 FALSE
## 3 2 7.685445 1.470564 3.198494 44 3.051178 3.419309 FALSE
##
## attr(,"class")
## [1] "gofOutlier"acg_dat["84", "Ins_conc_molal"] <- NA #replace outlier value with NA to exclude from analysis
plot_glx <- ggplot(acg_dat, aes(Glu_Gln_conc_molal))
plot_glx + geom_density()
plot_glx + geom_boxplot()
rosnerTest(acg_dat$Glu_Gln_conc_molal, k = 1)## $distribution
## [1] "Normal"
##
## $statistic
## R.1
## 2.840029
##
## $sample.size
## [1] 113
##
## $parameters
## k
## 1
##
## $alpha
## [1] 0.05
##
## $crit.value
## lambda.1
## 3.425263
##
## $n.outliers
## [1] 0
##
## $alternative
## [1] "Up to 1 observations are not\n from the same Distribution."
##
## $method
## [1] "Rosner's Test for Outliers"
##
## $data
## [1] 22.16208 24.54895 30.79902 25.15445 23.12144 23.45996 27.95219 28.85494
## [9] 28.47293 30.83024 24.02170 20.83895 26.81225 25.81936 25.83210 28.15351
## [17] 28.69984 24.45281 25.87441 26.10568 24.90431 27.69072 29.16664 25.73952
## [25] 31.00567 25.29233 29.35298 28.81394 30.45098 26.42377 25.25458 30.67838
## [33] 29.65833 34.25604 26.24353 28.06087 21.93668 25.00839 27.05607 29.92424
## [41] 27.12366 23.65348 26.82576 24.86136 26.92752 23.98651 26.21693 28.06806
## [49] 26.26637 27.44238 23.98526 25.40089 22.59948 23.46421 28.23542 33.70878
## [57] 26.99835 25.59893 26.95525 25.77541 32.34599 24.66776 25.46001 25.71687
## [65] 29.73266 25.54233 27.95401 29.18637 26.09421 27.74434 27.70571 21.88957
## [73] 26.72328 30.44184 21.32831 31.44932 28.73199 25.08766 29.57015 24.10930
## [81] 27.84161 24.85586 29.21009 20.93551 23.34771 27.38610 26.35670 26.10330
## [89] 28.90132 26.67329 24.94565 29.44010 29.32746 26.23233 24.32566 24.39736
## [97] 26.88103 30.91320 28.56651 27.47953 25.89977 24.22166 25.89081 28.12079
## [105] 29.77257 28.04112 29.21096 29.61668 27.08694 28.10033 24.05989 26.66364
## [113] 26.45077
##
## $data.name
## [1] "acg_dat$Glu_Gln_conc_molal"
##
## $bad.obs
## [1] 0
##
## $all.stats
## i Mean.i SD.i Value Obs.Num R.i+1 lambda.i+1 Outlier
## 1 0 26.82935 2.615006 34.25604 34 2.840029 3.425263 FALSE
##
## attr(,"class")
## [1] "gofOutlier"write.csv(acg_dat, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/acg_data_molal_outliers_removed_n113.csv", row.names = FALSE)LAG
plot_age <- ggplot(lag_dat, aes(mri_age_y))
plot_age + geom_density()
#age skewed
plot_naa <- ggplot(lag_dat, aes(NAA_conc_molal))
plot_naa + geom_density()
plot_naa + geom_boxplot()
rosnerTest(lag_dat$NAA_conc_molal, k = 10) #outliers val. 19.207, obs. 296 & Val. 10.148, obs. 306 - some issues in spectra, exclude from analysis for NAA## $distribution
## [1] "Normal"
##
## $statistic
## R.1 R.2 R.3 R.4 R.5 R.6 R.7 R.8
## 4.467031 4.044357 3.599759 3.613937 3.039885 2.720506 2.647000 2.651977
## R.9 R.10
## 2.660497 2.662447
##
## $sample.size
## [1] 321
##
## $parameters
## k
## 10
##
## $alpha
## [1] 0.05
##
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 lambda.5 lambda.6 lambda.7 lambda.8
## 3.742579 3.741706 3.740829 3.739949 3.739067 3.738181 3.737291 3.736399
## lambda.9 lambda.10
## 3.735503 3.734604
##
## $n.outliers
## [1] 2
##
## $alternative
## [1] "Up to 10 observations are not\n from the same Distribution."
##
## $method
## [1] "Rosner's Test for Outliers"
##
## $data
## [1] 14.16473 13.06076 14.35969 13.27563 16.92063 13.92045 14.90883 14.96705
## [9] 13.87538 13.79928 15.43119 13.74363 13.42774 14.32941 15.29816 15.36852
## [17] 14.96343 14.01462 15.11779 16.29892 15.28430 15.17137 15.34158 15.28575
## [25] 13.61378 16.76759 13.63702 14.41263 13.69061 15.64950 14.60942 15.25717
## [33] 15.68290 14.24180 15.96102 13.19084 14.22355 14.45688 13.89590 16.68967
## [41] 15.14115 13.88211 14.56337 14.77687 14.00542 14.12098 13.36743 13.70315
## [49] 14.59050 15.55276 13.33791 14.13365 12.08635 14.08027 13.90090 13.63249
## [57] 15.72428 11.86871 14.17563 11.96658 12.37485 14.57472 13.54554 14.35760
## [65] 13.86325 16.89273 14.68644 16.67785 15.93255 16.85334 14.93659 14.21829
## [73] 13.99673 14.23697 13.86957 13.86406 14.36716 14.41372 12.45632 14.76312
## [81] 14.75832 14.92673 13.88514 14.24089 15.42366 14.86462 14.32717 14.96667
## [89] 13.72426 14.16948 15.16894 13.93828 12.92563 13.94977 13.64160 14.55102
## [97] 14.51660 12.52392 14.21054 13.93670 16.19979 14.45536 14.10612 13.35969
## [105] 16.58862 14.41649 14.80811 14.26277 12.37256 15.44893 13.43728 13.57859
## [113] 13.72929 14.86513 11.95501 14.39891 14.11323 14.14168 14.73091 14.40429
## [121] 13.19741 13.71443 14.26528 14.92876 14.67063 14.35362 15.12298 13.95432
## [129] 13.68703 13.44733 14.01869 14.75568 14.74986 11.38458 12.08162 11.74260
## [137] 13.93366 14.16181 14.38636 14.40623 14.15278 14.17590 15.46643 14.09616
## [145] 13.62893 14.14878 13.54441 13.09883 12.79401 14.64345 14.02714 14.59131
## [153] 14.51512 15.41640 13.41648 13.54853 13.78080 13.79989 13.51420 14.77851
## [161] 15.60295 13.13390 13.99763 14.27699 15.00866 15.01634 12.09852 14.98452
## [169] 13.60870 16.34165 11.91384 12.33377 14.09992 15.06924 13.87378 13.47019
## [177] 14.02295 18.06089 14.25601 14.37746 13.24330 15.14502 13.95585 13.43714
## [185] 14.30145 16.30323 13.44700 14.70962 13.78765 13.43206 14.86922 14.73623
## [193] 14.54106 13.30455 14.61897 13.91250 15.44818 13.00756 14.21959 16.58017
## [201] 14.30707 15.14338 12.83493 13.96649 15.41109 13.77502 12.24289 15.19071
## [209] 13.32940 13.11120 15.00722 15.15742 13.96143 14.70099 14.27321 15.14983
## [217] 14.79365 15.19557 14.84745 13.19832 15.12751 13.95088 13.28456 14.90679
## [225] 12.94820 13.28282 14.39416 13.76430 14.91352 15.55618 15.20787 14.89958
## [233] 14.31622 14.53646 15.29391 13.51980 15.05685 14.73951 15.01094 15.61748
## [241] 14.09588 14.63147 15.59956 15.33848 14.07436 14.30913 13.94217 16.47833
## [249] 15.29907 13.85813 15.73137 15.08295 16.15169 14.28221 15.30824 14.47596
## [257] 14.47137 14.28477 15.33016 14.75400 14.54703 14.66072 13.49926 14.47973
## [265] 13.83620 14.83537 15.05792 14.18100 13.68113 13.28288 14.53224 13.92267
## [273] 14.80805 13.25470 14.60733 13.77756 15.36693 16.13917 14.36568 14.54767
## [281] 14.51644 16.31083 15.95003 13.41067 15.50889 13.79439 14.53462 14.83200
## [289] 14.42306 13.93431 13.79584 14.91428 12.66542 14.98849 15.38162 19.20747
## [297] 13.74790 17.99338 14.95991 13.93438 15.42154 13.27326 13.14415 14.47594
## [305] 13.99067 10.14812 13.73230 14.91913 15.85376 13.23125 13.40927 14.90803
## [313] 14.30314 14.54287 14.33013 16.55184 14.94880 14.30460 15.31501 14.46473
## [321] 16.31675
##
## $data.name
## [1] "lag_dat$NAA_conc_molal"
##
## $bad.obs
## [1] 0
##
## $all.stats
## i Mean.i SD.i Value Obs.Num R.i+1 lambda.i+1 Outlier
## 1 0 14.39149 1.0781156 19.20747 296 4.467031 3.742579 TRUE
## 2 1 14.37644 1.0454867 10.14812 306 4.044357 3.741706 TRUE
## 3 2 14.38970 1.0198436 18.06089 178 3.599759 3.740829 FALSE
## 4 3 14.37815 1.0003560 17.99338 298 3.613937 3.739949 FALSE
## 5 4 14.36675 0.9810139 11.38458 134 3.039885 3.739067 FALSE
## 6 5 14.37619 0.9680502 11.74260 136 2.720506 3.738181 FALSE
## 7 6 14.38455 0.9580956 16.92063 5 2.647000 3.737291 FALSE
## 8 7 14.37647 0.9488234 16.89273 66 2.651977 3.736399 FALSE
## 9 8 14.36843 0.9395707 11.86871 58 2.660497 3.735503 FALSE
## 10 9 14.37644 0.9303092 16.85334 70 2.662447 3.734604 FALSE
##
## attr(,"class")
## [1] "gofOutlier"lag_dat["296", "NAA_conc_molal"]<- NA
lag_dat["306", "NAA_conc_molal"]<- NA
plot_cre <- ggplot(lag_dat, aes(Cre_PCr_conc_molal))
plot_cre + geom_density()
plot_cre + geom_boxplot()
rosnerTest(lag_dat$Cre_PCr_conc_molal, k = 10) #outliers val. 14.829, obs. 241 (PS17_047) - structure in resid. & val. 14.09, obs. 178 (PS17_045, spectrum looked fine), excluding only PS17_047 for now## $distribution
## [1] "Normal"
##
## $statistic
## R.1 R.2 R.3 R.4 R.5 R.6 R.7 R.8
## 4.835902 4.220654 3.332561 3.278883 3.286432 3.307998 2.911341 2.895521
## R.9 R.10
## 2.767415 2.775535
##
## $sample.size
## [1] 321
##
## $parameters
## k
## 10
##
## $alpha
## [1] 0.05
##
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 lambda.5 lambda.6 lambda.7 lambda.8
## 3.742579 3.741706 3.740829 3.739949 3.739067 3.738181 3.737291 3.736399
## lambda.9 lambda.10
## 3.735503 3.734604
##
## $n.outliers
## [1] 2
##
## $alternative
## [1] "Up to 10 observations are not\n from the same Distribution."
##
## $method
## [1] "Rosner's Test for Outliers"
##
## $data
## [1] 9.979902 8.752796 9.874363 8.961297 11.510180 10.121568 10.418593
## [8] 11.028922 10.171413 9.929791 10.871830 11.064228 10.111309 9.319128
## [15] 10.662947 10.890644 10.197139 9.919287 10.426231 11.325982 10.557574
## [22] 11.088082 11.039646 11.172931 8.936929 12.109674 10.000368 10.420059
## [29] 9.572731 10.511121 9.984628 10.695378 11.039338 9.822119 11.097416
## [36] 9.239262 9.742404 10.192047 10.845538 12.117826 10.389433 9.535353
## [43] 10.010295 10.823639 9.876324 10.019161 8.608303 9.442680 11.444015
## [50] 11.419644 9.408142 10.139106 8.483950 9.792171 10.032519 8.968350
## [57] 11.562417 9.156420 10.285179 7.449145 10.159106 10.623228 10.245682
## [64] 11.514351 10.182653 11.786704 11.560836 12.034078 11.166365 11.463697
## [71] 10.397742 10.176308 11.796639 10.733740 10.174742 10.465407 9.824296
## [78] 11.420357 7.832529 10.284534 12.513826 11.576531 10.541650 10.808215
## [85] 10.483803 10.041962 10.382300 9.976776 10.765842 10.520082 11.422196
## [92] 10.843472 9.968723 10.781936 9.841352 9.930816 11.494602 9.908886
## [99] 10.765624 10.122441 11.914466 10.403873 10.056177 9.179081 11.252406
## [106] 10.011672 10.666794 9.981113 9.225683 10.367639 10.879950 11.335061
## [113] 9.425501 10.725975 10.241112 9.317008 9.889557 9.560765 10.593153
## [120] 11.263291 8.935236 10.137273 10.048855 11.211108 10.481657 10.028783
## [127] 11.460008 9.844337 10.207849 8.876822 9.557376 9.418484 9.631575
## [134] 11.004550 7.414051 8.984416 9.798072 10.236338 10.150536 10.108983
## [141] 10.653763 10.137920 10.541536 9.525275 9.886573 9.896683 9.185657
## [148] 9.642602 7.281217 9.968306 10.715037 10.739283 11.039381 11.250002
## [155] 9.759444 9.643949 9.675426 9.621214 9.141950 10.356422 10.788507
## [162] 9.362209 9.822128 9.645301 10.157630 10.381588 8.913720 11.109833
## [169] 9.890247 12.488886 9.115647 9.814255 10.627685 11.734993 10.295594
## [176] 10.173568 9.988731 14.090723 10.829617 10.247883 10.533434 10.999808
## [183] 10.001382 9.229915 10.354994 11.296388 9.994843 10.904378 9.345731
## [190] 9.618365 9.759703 10.530182 10.844392 8.895843 10.110129 11.423088
## [197] 11.096328 9.958072 9.931007 10.921273 9.154313 10.482202 10.423407
## [204] 10.629254 12.212044 9.913344 8.305599 10.446897 9.841379 9.379922
## [211] 9.844787 10.968766 10.801476 9.891839 10.632212 9.671280 10.872873
## [218] 10.666877 10.645958 9.493899 13.103099 9.512697 10.276565 9.846960
## [225] 9.773848 9.929345 11.034653 10.095542 10.422810 10.640722 10.407355
## [232] 11.450386 10.284795 10.179735 11.263605 10.183927 10.257883 10.937157
## [239] 10.529431 11.012918 14.829071 9.133491 11.531961 10.173887 9.610385
## [246] 9.349643 9.324653 10.926724 11.044735 8.859819 11.228731 10.816240
## [253] 11.271292 9.861827 10.852278 10.764773 9.380198 9.481786 10.781017
## [260] 8.985104 9.888393 9.434081 8.870117 9.122209 9.542512 9.889146
## [267] 10.769356 9.616033 8.966342 9.255494 10.764599 9.897261 9.427714
## [274] 8.163350 10.017130 9.265501 10.944724 12.004282 9.764853 9.856451
## [281] 10.020483 12.009948 10.044769 9.667192 11.197137 9.964094 10.035066
## [288] 10.812708 9.728669 9.918824 8.843802 10.202045 9.183452 9.822267
## [295] 10.790509 12.454846 8.641534 9.507122 10.361671 9.646740 11.011202
## [302] 9.075410 10.166484 9.737467 9.615290 7.882414 10.377305 10.224366
## [309] 12.449600 10.083783 9.260131 11.009901 9.914846 9.659798 10.610323
## [316] 11.561746 11.077058 9.291690 10.160005 9.934510 11.934339
##
## $data.name
## [1] "lag_dat$Cre_PCr_conc_molal"
##
## $bad.obs
## [1] 0
##
## $all.stats
## i Mean.i SD.i Value Obs.Num R.i+1 lambda.i+1 Outlier
## 1 0 10.26105 0.9446048 14.829071 241 4.835902 3.742579 TRUE
## 2 1 10.24678 0.9107457 14.090723 178 4.220654 3.741706 TRUE
## 3 2 10.23473 0.8862592 7.281217 149 3.332561 3.740829 FALSE
## 4 3 10.24402 0.8719682 13.103099 221 3.278883 3.739949 FALSE
## 5 4 10.23500 0.8583617 7.414051 135 3.286432 3.739067 FALSE
## 6 5 10.24393 0.8448557 7.449145 60 3.307998 3.738181 FALSE
## 7 6 10.25280 0.8313244 7.832529 79 2.911341 3.737291 FALSE
## 8 7 10.26051 0.8213002 7.882414 306 2.895521 3.736399 FALSE
## 9 8 10.26810 0.8114875 12.513826 81 2.767415 3.735503 FALSE
## 10 9 10.26091 0.8027211 12.488886 170 2.775535 3.734604 FALSE
##
## attr(,"class")
## [1] "gofOutlier"lag_dat["241", "Cre_PCr_conc_molal"]<-NA #replace outlier value with NA to exclude from analysis
plot_cho <- ggplot(lag_dat, aes(Cho_GPC_PCh_conc_molal))
plot_cho + geom_density()
plot_cho + geom_boxplot()
rosnerTest(lag_dat$Cho_GPC_PCh_conc_molal, k = 8)## $distribution
## [1] "Normal"
##
## $statistic
## R.1 R.2 R.3 R.4 R.5 R.6 R.7 R.8
## 3.694967 3.121637 3.061660 3.062209 2.979166 3.001536 2.871646 2.733239
##
## $sample.size
## [1] 321
##
## $parameters
## k
## 8
##
## $alpha
## [1] 0.05
##
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 lambda.5 lambda.6 lambda.7 lambda.8
## 3.742579 3.741706 3.740829 3.739949 3.739067 3.738181 3.737291 3.736399
##
## $n.outliers
## [1] 0
##
## $alternative
## [1] "Up to 8 observations are not\n from the same Distribution."
##
## $method
## [1] "Rosner's Test for Outliers"
##
## $data
## [1] 2.068181 2.062467 2.386391 2.355035 2.498191 2.256110 2.077815 2.425292
## [9] 2.398770 1.895982 1.865260 2.184158 2.004899 2.128662 2.236830 2.378977
## [17] 2.355458 2.266506 2.404117 2.307009 2.052566 2.364941 2.267604 1.940308
## [25] 2.587651 2.716570 2.469831 2.547177 2.087747 2.170207 2.232363 2.241430
## [33] 2.441839 2.264953 2.333662 2.359801 2.270399 2.241648 1.920730 2.467142
## [41] 2.242001 2.227010 2.229581 2.053798 1.939191 2.364334 2.678152 2.362113
## [49] 1.676867 2.166236 2.306276 2.211050 2.178445 2.361600 2.344510 2.268087
## [57] 2.667157 2.059386 2.034573 2.045022 2.214431 2.995506 1.846508 2.183946
## [65] 2.072392 2.429630 2.330074 2.491408 2.265040 2.154392 2.411001 2.046897
## [73] 2.534412 2.027012 2.376839 2.222266 2.306644 2.292419 2.108096 2.172672
## [81] 1.839589 1.859195 1.443755 1.640698 2.420001 2.337525 2.238790 2.240345
## [89] 2.290795 2.399979 2.350376 2.151479 2.142910 2.110996 1.985190 1.928934
## [97] 2.304258 1.765061 2.137595 1.985591 2.781876 2.123475 1.756957 1.700318
## [105] 1.975316 1.974756 1.999963 2.031005 1.687589 1.708882 2.296279 2.182991
## [113] 2.276755 2.101586 1.837847 2.371711 2.065515 2.499424 2.363551 2.223515
## [121] 2.232809 2.045991 1.947255 2.339091 2.193248 2.528316 2.691920 2.777161
## [129] 2.451074 2.763753 2.377191 2.160428 2.215103 2.214312 2.533792 2.384939
## [137] 2.179386 2.552211 2.673232 2.476304 1.625876 1.926042 2.202001 2.204096
## [145] 2.344069 2.152159 1.895155 1.710908 1.401447 1.946252 2.047116 2.456831
## [153] 2.086612 2.201666 2.190181 2.294616 2.581968 2.225911 2.498999 2.000618
## [161] 2.163682 1.960245 2.204252 2.159811 2.097381 2.254865 1.844443 1.598729
## [169] 1.951488 3.208685 1.836756 1.950559 1.992121 2.281328 2.008449 2.599790
## [177] 2.317712 3.033330 2.698335 2.332814 2.583764 2.234143 2.191218 2.264957
## [185] 1.977277 2.067888 2.356329 2.431505 2.779803 2.150897 2.566469 2.101205
## [193] 2.159059 2.164490 2.152499 2.489624 2.791537 1.673365 2.193421 2.761471
## [201] 1.959300 2.022674 1.970196 2.366537 2.774857 2.129388 1.836375 2.360163
## [209] 2.476086 2.665161 2.258604 2.298965 2.012883 2.201029 2.158646 2.353116
## [217] 1.946114 1.912304 1.692442 2.307008 2.449170 2.132806 2.045527 2.155850
## [225] 2.166616 2.440956 1.783126 2.272129 2.595674 2.284913 2.155251 2.400964
## [233] 2.048543 2.045160 2.356748 2.222524 2.346513 2.378440 2.281038 1.906127
## [241] 1.725583 2.479767 2.262261 2.082256 2.047599 2.138599 2.252565 2.262486
## [249] 1.881557 2.322610 1.489680 2.225178 2.118903 2.390375 2.437536 2.400945
## [257] 1.832693 2.290250 2.007843 2.346441 2.089229 2.428740 2.391081 2.732220
## [265] 1.971298 2.235255 2.116196 1.941129 2.009743 1.999019 1.927678 1.994412
## [273] 2.372282 2.491822 1.886355 2.109659 2.028203 1.807961 2.439615 2.152157
## [281] 1.892138 2.163111 2.376896 1.914835 2.734225 1.944843 1.756721 2.135482
## [289] 2.101386 1.768588 1.974516 1.975407 2.167757 2.112266 2.510332 2.521732
## [297] 2.431779 2.959724 2.046249 1.986009 2.300110 2.368720 2.151772 2.070465
## [305] 2.457353 2.243696 2.105342 2.108232 2.180180 2.349512 1.911659 2.150353
## [313] 2.011103 2.243254 2.417149 1.800195 2.190437 2.875345 2.116343 2.319437
## [321] 2.422916
##
## $data.name
## [1] "lag_dat$Cho_GPC_PCh_conc_molal"
##
## $bad.obs
## [1] 0
##
## $all.stats
## i Mean.i SD.i Value Obs.Num R.i+1 lambda.i+1 Outlier
## 1 0 2.208213 0.2707661 3.208685 170 3.694967 3.742579 FALSE
## 2 1 2.205087 0.2653235 3.033330 178 3.121637 3.741706 FALSE
## 3 2 2.202490 0.2616371 1.401447 149 3.061660 3.740829 FALSE
## 4 3 2.205009 0.2581460 2.995506 62 3.062209 3.739949 FALSE
## 5 4 2.202516 0.2546890 1.443755 83 2.979166 3.739067 FALSE
## 6 5 2.204917 0.2514736 2.959724 298 3.001536 3.738181 FALSE
## 7 6 2.202521 0.2482341 1.489680 251 2.871646 3.737291 FALSE
## 8 7 2.204791 0.2453333 2.875345 318 2.733239 3.736399 FALSE
##
## attr(,"class")
## [1] "gofOutlier"plot_ins <- ggplot(lag_dat, aes(Ins_conc_molal))
plot_ins + geom_density()
plot_ins + geom_boxplot()
rosnerTest(lag_dat$Ins_conc_molal, k = 8) #outliers: val. 14.077, obs. 241 (PS17_047) - structure in resid; val. 11.777, obs. 221 (PS16_014) - spectrum fine, keep for now; val. 1.542, obs. 197 (PS18_028) - noisy spectrum, spiky residuals, exclude; val. 2.04, obs. 62 (PS14_111) - spectrum fine, keep for now## $distribution
## [1] "Normal"
##
## $statistic
## R.1 R.2 R.3 R.4 R.5 R.6 R.7 R.8
## 5.981441 4.463363 4.077869 3.761599 3.384729 3.303295 3.148162 3.105675
##
## $sample.size
## [1] 321
##
## $parameters
## k
## 8
##
## $alpha
## [1] 0.05
##
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 lambda.5 lambda.6 lambda.7 lambda.8
## 3.742579 3.741706 3.740829 3.739949 3.739067 3.738181 3.737291 3.736399
##
## $n.outliers
## [1] 4
##
## $alternative
## [1] "Up to 8 observations are not\n from the same Distribution."
##
## $method
## [1] "Rosner's Test for Outliers"
##
## $data
## [1] 8.084411 7.069186 6.913031 6.915494 7.726417 8.046174 6.937152
## [8] 4.793633 5.952276 6.984233 6.400936 4.413732 7.176021 7.448535
## [15] 7.547442 6.541306 8.400250 5.957678 6.965664 7.688786 6.245160
## [22] 7.360921 7.925104 4.710055 6.035369 7.937845 7.166536 7.269577
## [29] 5.736145 7.751891 6.217281 7.121735 7.138640 6.104017 7.736566
## [36] 6.539681 4.958196 7.275083 6.485777 6.598573 7.153802 7.352086
## [43] 7.119035 6.146843 6.038556 6.638197 4.629802 6.187281 6.553798
## [50] 5.945866 7.474194 5.559932 5.743412 6.163849 8.347824 4.432658
## [57] 8.169688 5.045280 8.873237 6.213139 6.400722 2.040275 8.158702
## [64] 6.760986 6.468052 7.273556 6.742605 7.552488 7.843304 5.876394
## [71] 8.162408 4.963944 7.945325 5.440451 4.854219 4.834565 5.839066
## [78] 5.839637 5.890621 5.871859 6.963936 5.579543 5.190653 7.062975
## [85] 6.065434 7.046003 6.865074 6.210505 6.024185 6.121049 8.133551
## [92] 6.863480 6.537634 7.457231 5.008241 6.892333 5.936366 3.826364
## [99] 7.075620 5.567096 5.248858 5.491554 4.529573 7.011630 5.896171
## [106] 6.268198 5.292489 5.663161 6.753982 5.643647 6.788746 5.437125
## [113] 3.078293 7.511207 2.568354 8.286250 5.309958 6.053735 8.092627
## [120] 7.122038 6.509101 6.564424 5.780997 6.770400 5.727577 7.829891
## [127] 6.206182 7.258719 6.558572 4.787960 6.132452 6.029208 6.663446
## [134] 7.630692 4.608510 5.731585 5.435869 5.209567 2.734199 5.857696
## [141] 5.638088 6.135153 6.888234 5.236533 5.655180 4.948126 4.992038
## [148] 6.025793 2.972647 6.870427 7.728083 7.800580 7.334580 6.676194
## [155] 6.170230 5.196197 7.073851 6.127912 7.232895 4.960469 6.612268
## [162] 6.088183 5.798318 6.543746 6.997835 6.987522 4.250557 4.637685
## [169] 6.778603 7.502499 4.679992 5.270903 6.367636 6.301188 7.918660
## [176] 6.533953 7.497948 8.801777 6.797100 5.264895 5.014240 6.141243
## [183] 5.472028 3.885218 4.816043 6.248540 4.626087 6.007857 6.538200
## [190] 7.586651 7.411192 7.320786 7.682231 7.120925 6.686811 7.373729
## [197] 1.541994 6.314608 6.840831 7.095049 6.617076 6.056406 5.882616
## [204] 4.775951 5.731427 5.979244 6.418407 6.735628 5.243289 5.869558
## [211] 5.226573 6.275003 7.093997 7.541310 5.731832 4.649731 6.150065
## [218] 6.428658 7.426624 6.729085 11.777066 5.591387 3.080430 8.027587
## [225] 7.066068 6.661947 6.600256 6.955873 6.949031 6.824263 7.605754
## [232] 6.414457 6.094616 6.805794 5.865261 8.270934 4.935241 6.818562
## [239] 5.441384 7.674150 14.076805 4.939637 6.137708 5.300423 6.306576
## [246] 6.245176 7.663984 7.103391 8.335725 3.649002 9.145931 7.122485
## [253] 8.331023 7.199990 8.067478 7.557525 5.521439 6.036198 6.082461
## [260] 7.524328 5.926927 5.705055 5.399996 6.232379 5.707127 4.605042
## [267] 4.333050 5.607343 5.940738 5.740379 6.177392 6.685727 5.210419
## [274] 6.865329 6.749236 5.464544 6.496419 6.705573 7.274507 6.221882
## [281] 5.994352 6.059500 7.660503 6.660801 7.439851 6.777698 5.282780
## [288] 8.773125 6.173129 5.673978 5.639679 5.174471 8.615918 5.774767
## [295] 5.411721 7.832813 7.925781 4.996084 6.276000 6.161567 6.484734
## [302] 5.594141 6.793109 6.401047 9.077790 8.562826 5.674753 6.368745
## [309] 8.425749 6.658481 4.433254 5.482190 7.149805 6.079902 6.139307
## [316] 7.111624 5.971685 6.120964 6.571466 5.418973 6.379001
##
## $data.name
## [1] "lag_dat$Ins_conc_molal"
##
## $bad.obs
## [1] 0
##
## $all.stats
## i Mean.i SD.i Value Obs.Num R.i+1 lambda.i+1 Outlier
## 1 0 6.383077 1.286267 14.076805 241 5.981441 3.742579 TRUE
## 2 1 6.359034 1.213890 11.777066 221 4.463363 3.741706 TRUE
## 3 2 6.342050 1.177099 1.541994 197 4.077869 3.740829 TRUE
## 4 3 6.357144 1.147615 2.040275 62 3.761599 3.739949 TRUE
## 5 4 6.370762 1.123401 2.568354 115 3.384729 3.739067 FALSE
## 6 5 6.382795 1.104533 2.734199 139 3.303295 3.738181 FALSE
## 7 6 6.394378 1.086898 2.972647 149 3.148162 3.737291 FALSE
## 8 7 6.405275 1.071259 3.078293 113 3.105675 3.736399 FALSE
##
## attr(,"class")
## [1] "gofOutlier"lag_dat["241", "Ins_conc_molal"] <- NA #replace outlier value with NA to exclude from analysis
lag_dat["197", "Ins_conc_molal"] <- NA
plot_glx <- ggplot(lag_dat, aes(Glu_Gln_conc_molal))
plot_glx + geom_density()
plot_glx + geom_boxplot()
rosnerTest(lag_dat$Glu_Gln_conc_molal, k = 7)#outliers: val. 10.09, obs. 241 (PS17_047) - structure in resid; val. 31.2899, obs. 296 (PS21_030) - big spike in resid. around 2.1, exclude.## $distribution
## [1] "Normal"
##
## $statistic
## R.1 R.2 R.3 R.4 R.5 R.6 R.7
## 4.954047 4.247538 3.396313 3.439741 3.463175 3.267721 3.232088
##
## $sample.size
## [1] 321
##
## $parameters
## k
## 7
##
## $alpha
## [1] 0.05
##
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 lambda.5 lambda.6 lambda.7
## 3.742579 3.741706 3.740829 3.739949 3.739067 3.738181 3.737291
##
## $n.outliers
## [1] 2
##
## $alternative
## [1] "Up to 7 observations are not\n from the same Distribution."
##
## $method
## [1] "Rosner's Test for Outliers"
##
## $data
## [1] 23.58311 20.72833 21.43465 20.79838 24.55542 19.38008 22.13533 21.98148
## [9] 18.58332 21.80094 23.39061 22.39400 22.25994 19.84839 22.62240 19.74196
## [17] 19.55718 20.18644 20.93614 24.01329 22.14911 23.96545 23.49753 23.68093
## [25] 20.83139 25.55075 18.83744 23.02848 20.10324 24.84118 23.11824 25.62522
## [33] 25.14532 20.60338 23.10183 18.68727 19.90887 22.63726 23.44172 23.97287
## [41] 20.79684 21.03860 21.41871 22.38334 19.07457 19.80466 21.54623 21.54133
## [49] 22.80842 25.04145 20.65555 22.02043 18.52301 21.98565 22.88652 20.99602
## [57] 23.64194 20.45706 24.44658 22.40116 19.22839 24.77315 20.40702 22.32472
## [65] 22.79321 23.27919 23.66026 23.59194 23.28915 23.15967 20.72270 22.14974
## [73] 23.87191 27.09158 23.12072 25.11318 22.08527 21.17324 14.25503 21.44083
## [81] 23.08444 21.47161 23.79592 26.39922 24.52352 23.18206 22.57951 22.53444
## [89] 24.24539 23.38651 21.39081 23.50524 21.84539 23.16617 20.66627 20.68771
## [97] 21.33913 23.82111 24.17549 20.43359 24.17777 21.80255 21.79162 21.34246
## [105] 26.43216 21.02534 21.61161 22.57479 21.58113 22.47585 21.89183 22.86451
## [113] 21.63056 23.83333 23.16833 22.80824 21.74170 19.79860 22.68433 25.43932
## [121] 20.46600 19.55287 21.71921 23.91575 19.94617 22.54020 25.12373 20.26601
## [129] 20.42848 18.40115 22.31740 21.34145 21.18639 24.46723 17.55526 19.17975
## [137] 22.14220 21.22080 20.62138 19.29726 22.33080 22.36491 23.64798 24.04583
## [145] 25.02751 22.02316 20.84404 22.82850 23.75262 19.93399 21.93649 21.46368
## [153] 23.51970 23.79151 22.29261 22.55887 19.88319 23.52534 18.26039 22.44651
## [161] 23.87196 21.33171 21.19221 21.21896 25.13735 21.72890 21.00404 24.38403
## [169] 22.77448 24.07955 20.92792 22.79597 21.53526 25.58587 20.23221 20.38536
## [177] 19.64011 29.12643 21.85437 22.81675 24.27197 20.09271 21.28381 21.20096
## [185] 20.61092 26.55674 20.24095 22.22341 21.17598 21.05662 20.70603 22.17497
## [193] 22.22738 17.77073 21.62541 23.02400 23.16615 17.38515 20.64300 21.97048
## [201] 23.41849 22.42124 20.67899 21.35918 23.92414 18.75710 19.72431 22.35249
## [209] 23.13971 19.66645 21.97753 21.86935 22.15373 18.04793 20.63749 19.81250
## [217] 22.91987 21.99328 21.86022 19.23943 22.18618 21.97414 22.10059 21.96999
## [225] 21.38803 21.08812 21.12026 20.43893 20.01495 22.48505 23.44854 23.27656
## [233] 22.72471 21.03119 22.44160 19.23947 24.24070 21.72010 21.19354 21.24039
## [241] 10.09497 19.32810 23.68178 22.11865 18.91616 19.59606 19.02554 22.92979
## [249] 26.99382 18.91880 22.36158 21.71252 22.51005 20.90044 23.69281 22.60270
## [257] 20.69658 19.00693 19.97580 15.08181 17.70438 17.08767 17.40514 16.78328
## [265] 20.40818 19.57305 21.57031 21.73287 18.10301 16.64465 21.45075 20.54401
## [273] 19.32920 14.36555 21.66594 17.20659 21.57909 24.57946 17.49328 18.35123
## [281] 21.22156 24.11680 22.51873 19.19853 23.08553 19.03798 23.38476 21.59376
## [289] 19.09657 21.28636 18.16122 19.70762 18.41625 18.73419 22.65583 31.28994
## [297] 18.89292 17.32671 18.83142 19.60466 24.66522 18.35557 22.92964 18.49429
## [305] 21.56271 18.98307 20.17174 24.36413 25.12361 21.76663 19.89145 21.01574
## [313] 20.06260 25.41609 26.15708 23.69646 23.29001 20.54746 23.65148 20.44966
## [321] 28.51150
##
## $data.name
## [1] "lag_dat$Glu_Gln_conc_molal"
##
## $bad.obs
## [1] 0
##
## $all.stats
## i Mean.i SD.i Value Obs.Num R.i+1 lambda.i+1 Outlier
## 1 0 21.68857 2.340229 10.09497 241 4.954047 3.742579 TRUE
## 2 1 21.72480 2.251926 31.28994 296 4.247538 3.741706 TRUE
## 3 2 21.69482 2.190549 14.25503 79 3.396313 3.740829 FALSE
## 4 3 21.71821 2.153715 29.12643 178 3.439741 3.739949 FALSE
## 5 4 21.69484 2.116351 14.36555 274 3.463175 3.739067 FALSE
## 6 5 21.71804 2.078963 28.51150 321 3.267721 3.738181 FALSE
## 7 6 21.69647 2.046559 15.08181 260 3.232088 3.737291 FALSE
##
## attr(,"class")
## [1] "gofOutlier"lag_dat["241", "Glu_Gln_conc_molal"] <- NA
lag_dat["296", "Glu_Gln_conc_molal"] <- NA
#exclude PS17_047 from all analyses due to outlier status on multiple metabolites and issues in spectrum (noisy, structure in resid)
lag_dat["241", "NAA_conc_molal"] <- NA
lag_dat["241", "Cre_PCr_conc_molal"] <- NA
lag_dat["241", "Cho_GPC_PCh_conc_molal"] <- NA
write.csv(lag_dat, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/lag_data_molal_outliers_removed_n320.csv", row.names = FALSE)##Descriptive stats
#ACG - count individual subject N and N of females
acg_sub_sex <- acg_dat %>%
dplyr::select(subj_id, female)
acg_unique_subs<-unique(acg_sub_sex)
length(acg_unique_subs$subj_id)## [1] 102sum(acg_unique_subs$female)## [1] 47describe(acg_dat)## vars n mean sd median trimmed mad min
## study_code* 1 113 57.00 32.76 57.00 57.00 41.51 1.00
## subj_id* 2 113 51.42 28.89 54.00 51.46 35.58 1.00
## spec_id 3 113 5072.42 3034.48 3244.00 4708.09 1254.28 2058.00
## date_scan* 4 113 49.21 28.71 50.00 49.32 37.06 1.00
## hand* 5 113 4.54 1.04 5.00 4.78 0.00 1.00
## female 6 113 0.47 0.50 0.00 0.46 0.00 0.00
## mri_age_y 7 113 4.06 1.10 3.68 3.97 1.00 2.34
## folder* 8 113 57.00 32.76 57.00 57.00 41.51 1.00
## File_ID* 9 113 57.00 32.76 57.00 57.00 41.51 1.00
## fGM 10 113 0.77 0.05 0.78 0.78 0.05 0.55
## fWM 11 113 0.04 0.03 0.03 0.03 0.02 0.00
## fCSF 12 113 0.19 0.05 0.18 0.19 0.04 0.09
## NAA_conc_molal 13 113 14.89 1.32 15.03 14.94 1.25 11.24
## Cre_PCr_conc_molal 14 113 12.04 1.07 12.04 12.06 0.84 8.41
## Cho_GPC_PCh_conc_molal 15 112 2.50 0.34 2.48 2.49 0.30 1.57
## Ins_conc_molal 16 112 7.64 1.53 7.71 7.72 1.46 2.99
## Glu_conc_molal 17 113 22.38 2.25 22.50 22.46 2.01 16.03
## Glu_Gln_conc_molal 18 113 26.83 2.62 26.72 26.81 2.64 20.84
## max range skew kurtosis se
## study_code* 113.00 112.00 0.00 -1.23 3.08
## subj_id* 102.00 101.00 -0.05 -1.15 2.72
## spec_id 11448.00 9390.00 0.87 -0.81 285.46
## date_scan* 97.00 96.00 -0.04 -1.29 2.70
## hand* 6.00 5.00 -2.07 3.38 0.10
## female 1.00 1.00 0.12 -2.00 0.05
## mri_age_y 7.36 5.03 0.72 -0.49 0.10
## folder* 113.00 112.00 0.00 -1.23 3.08
## File_ID* 113.00 112.00 0.00 -1.23 3.08
## fGM 0.87 0.32 -1.14 1.59 0.01
## fWM 0.29 0.29 4.07 25.28 0.00
## fCSF 0.36 0.27 0.90 0.92 0.00
## NAA_conc_molal 18.30 7.05 -0.37 0.21 0.12
## Cre_PCr_conc_molal 14.76 6.35 -0.36 0.83 0.10
## Cho_GPC_PCh_conc_molal 3.33 1.76 0.16 0.08 0.03
## Ins_conc_molal 10.89 7.90 -0.51 0.35 0.14
## Glu_conc_molal 28.20 12.17 -0.31 0.40 0.21
## Glu_Gln_conc_molal 34.26 13.42 0.13 -0.01 0.25#LAG - count individual subject N and N of females
lag_sub_sex <- lag_dat %>%
dplyr::select(subj_id, female)
lag_unique_subs<-unique(lag_sub_sex)
length(lag_unique_subs$subj_id)## [1] 95sum(lag_unique_subs$female)## [1] 47describe(lag_dat)## vars n mean sd median trimmed mad min
## study_code* 1 321 160.94 92.72 161.00 161.00 118.61 1.00
## subj_id* 2 321 40.27 24.55 41.00 39.16 28.17 1.00
## spec_id 3 320 7744.26 3536.93 6379.50 7565.04 3846.61 3118.00
## date_scan* 4 321 143.73 82.33 143.00 143.92 103.78 1.00
## hand* 5 321 4.64 0.90 5.00 4.88 0.00 1.00
## female 6 321 0.47 0.50 0.00 0.46 0.00 0.00
## mri_age_y 7 321 5.95 1.90 5.39 5.79 1.79 2.41
## folder* 8 321 161.00 92.81 161.00 161.00 118.61 1.00
## File_ID* 9 321 161.00 92.81 161.00 161.00 118.61 1.00
## fGM 10 321 0.62 0.10 0.64 0.63 0.10 0.30
## fWM 11 321 0.35 0.11 0.34 0.34 0.11 0.00
## fCSF 12 321 0.03 0.02 0.03 0.03 0.01 0.00
## NAA_conc_molal 13 318 14.39 1.02 14.36 14.37 0.88 11.38
## Cre_PCr_conc_molal 14 320 10.25 0.91 10.17 10.23 0.83 7.28
## Cho_GPC_PCh_conc_molal 15 320 2.21 0.27 2.20 2.20 0.24 1.40
## Ins_conc_molal 16 319 6.37 1.19 6.37 6.39 1.08 2.04
## Glu_conc_molal 17 321 17.88 1.90 17.87 17.88 1.73 10.63
## Glu_Gln_conc_molal 18 319 21.69 2.19 21.74 21.72 2.02 14.26
## max range skew kurtosis se
## study_code* 320.00 319.00 0.00 -1.21 5.17
## subj_id* 95.00 94.00 0.28 -0.83 1.37
## spec_id 14205.00 11087.00 0.35 -1.37 197.72
## date_scan* 284.00 283.00 0.00 -1.19 4.60
## hand* 6.00 5.00 -2.35 4.48 0.05
## female 1.00 1.00 0.13 -1.99 0.03
## mri_age_y 11.13 8.72 0.66 -0.46 0.11
## folder* 321.00 320.00 0.00 -1.21 5.18
## File_ID* 321.00 320.00 0.00 -1.21 5.18
## fGM 0.80 0.50 -0.75 0.40 0.01
## fWM 0.66 0.65 0.43 0.20 0.01
## fCSF 0.21 0.21 3.12 16.87 0.00
## NAA_conc_molal 18.06 6.68 0.20 0.95 0.06
## Cre_PCr_conc_molal 14.09 6.81 0.17 1.33 0.05
## Cho_GPC_PCh_conc_molal 3.21 1.81 0.22 0.80 0.02
## Ins_conc_molal 11.78 9.74 -0.13 1.69 0.07
## Glu_conc_molal 24.87 14.24 -0.07 1.22 0.11
## Glu_Gln_conc_molal 29.13 14.87 -0.12 0.74 0.12#Calculate total Ns from both datasets combined
all_dat<-bind_rows(acg_dat, lag_dat)
length(unique(all_dat$File_ID))## [1] 434length(unique(all_dat$subj_id))## [1] 127all_sub_sex <- all_dat %>%
dplyr::select(subj_id, female)
all_unique_subs<-unique(all_sub_sex)
sum(all_unique_subs$female)## [1] 62describe(all_dat)## vars n mean sd median trimmed mad min
## study_code* 1 434 216.94 124.93 216.50 216.97 160.12 1.00
## subj_id* 2 434 53.87 33.15 54.00 52.20 38.55 1.00
## spec_id 3 433 7046.99 3606.23 5832.00 6818.88 3964.47 2058.00
## date_scan* 4 434 184.27 105.60 183.00 184.45 133.43 1.00
## hand* 5 434 4.62 0.93 5.00 4.86 0.00 1.00
## female 6 434 0.47 0.50 0.00 0.46 0.00 0.00
## mri_age_y 7 434 5.46 1.92 4.95 5.28 1.68 2.34
## folder* 8 434 217.04 125.00 217.50 217.05 160.12 1.00
## File_ID* 9 434 217.50 125.43 217.50 217.50 160.86 1.00
## fGM 10 434 0.66 0.11 0.67 0.67 0.11 0.30
## fWM 11 434 0.27 0.17 0.29 0.27 0.16 0.00
## fCSF 12 434 0.07 0.08 0.03 0.06 0.03 0.00
## NAA_conc_molal 13 431 14.52 1.13 14.48 14.51 0.99 11.24
## Cre_PCr_conc_molal 14 433 10.71 1.24 10.54 10.65 1.10 7.28
## Cho_GPC_PCh_conc_molal 15 432 2.28 0.32 2.26 2.27 0.26 1.40
## Ins_conc_molal 16 431 6.70 1.40 6.66 6.68 1.28 2.04
## Glu_conc_molal 17 434 19.05 2.81 18.53 18.89 2.53 10.63
## Glu_Gln_conc_molal 18 432 23.04 3.23 22.54 22.85 2.77 14.26
## max range skew kurtosis se
## study_code* 432.00 431.00 0.00 -1.21 6.00
## subj_id* 127.00 126.00 0.33 -0.78 1.59
## spec_id 14205.00 12147.00 0.45 -1.19 173.30
## date_scan* 366.00 365.00 0.00 -1.18 5.07
## hand* 6.00 5.00 -2.29 4.27 0.04
## female 1.00 1.00 0.13 -1.99 0.02
## mri_age_y 11.13 8.79 0.79 -0.09 0.09
## folder* 433.00 432.00 0.00 -1.20 6.00
## File_ID* 434.00 433.00 0.00 -1.21 6.02
## fGM 0.87 0.57 -0.60 0.10 0.01
## fWM 0.66 0.66 -0.09 -0.87 0.01
## fCSF 0.36 0.36 1.37 0.83 0.00
## NAA_conc_molal 18.30 7.05 0.10 0.57 0.05
## Cre_PCr_conc_molal 14.76 7.48 0.44 0.20 0.06
## Cho_GPC_PCh_conc_molal 3.33 1.93 0.45 0.70 0.02
## Ins_conc_molal 11.78 9.74 0.08 0.70 0.07
## Glu_conc_molal 28.20 17.57 0.49 0.13 0.13
## Glu_Gln_conc_molal 34.26 20.00 0.54 0.34 0.16LME modelling ACG
run mixed-effects models for ACG with participant as random effect, age, sex and age x sex interaction as fixed effects NAA run on acg_data_molal_outliers_removed_n113.csv
acg_naa_age<- lmer(NAA_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_naa_age)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## Data: acg_dat
##
## AIC BIC logLik deviance df.resid
## 383.3 394.2 -187.6 375.3 109
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.53495 -0.39238 0.05273 0.38044 1.68117
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 1.1187 1.0577
## Residual 0.6121 0.7824
## Number of obs: 113, groups: subj_id, 102
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 13.8379 0.4588 108.1511 30.16 <2e-16 ***
## mri_age_y 0.2513 0.1079 105.1446 2.33 0.0217 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## mri_age_y -0.959acg_naa_age_sex<- lmer(NAA_conc_molal ~ mri_age_y + female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_naa_age_sex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + female + (1 | subj_id)
## Data: acg_dat
##
## AIC BIC logLik deviance df.resid
## 385.2 398.8 -187.6 375.2 108
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.54531 -0.40334 0.05382 0.37656 1.67013
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 1.1189 1.0578
## Residual 0.6113 0.7819
## Number of obs: 113, groups: subj_id, 102
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 13.85437 0.46506 109.75797 29.791 <2e-16 ***
## mri_age_y 0.25366 0.10837 104.64089 2.341 0.0211 *
## female -0.05648 0.26075 99.22276 -0.217 0.8289
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y
## mri_age_y -0.926
## female -0.165 -0.098anova(acg_naa_age,acg_naa_age_sex)## Data: acg_dat
## Models:
## acg_naa_age: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_naa_age_sex: NAA_conc_molal ~ mri_age_y + female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## acg_naa_age 4 383.25 394.16 -187.62 375.25
## acg_naa_age_sex 5 385.20 398.84 -187.60 375.20 0.0469 1 0.8285acg_naa_ageXsex<- lmer(NAA_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_naa_ageXsex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 |
## subj_id)
## Data: acg_dat
##
## AIC BIC logLik deviance df.resid
## 387.1 403.4 -187.5 375.1 107
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.55254 -0.39337 0.03275 0.37793 1.69016
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 1.1179 1.0573
## Residual 0.6105 0.7813
## Number of obs: 113, groups: subj_id, 102
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 13.6886 0.6657 106.4622 20.562 <2e-16 ***
## mri_age_y 0.2954 0.1616 107.8036 1.827 0.0704 .
## female 0.2522 0.9254 110.1751 0.273 0.7857
## mri_age_y:female -0.0757 0.2178 109.0452 -0.348 0.7288
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y female
## mri_age_y -0.964
## female -0.719 0.694
## mri_g_y:fml 0.716 -0.742 -0.960anova(acg_naa_age,acg_naa_ageXsex)## Data: acg_dat
## Models:
## acg_naa_age: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_naa_ageXsex: NAA_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## acg_naa_age 4 383.25 394.16 -187.62 375.25
## acg_naa_ageXsex 6 387.08 403.45 -187.54 375.08 0.1677 2 0.9196#including sex in model does not significantly improve fit Cre
acg_cre_age<- lmer(Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_cre_age)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
## Data: acg_dat
##
## AIC BIC logLik deviance df.resid
## 340.9 351.8 -166.5 332.9 109
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.12075 -0.39494 0.02556 0.37255 1.62021
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.6867 0.8287
## Residual 0.4877 0.6983
## Number of obs: 113, groups: subj_id, 102
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 11.67366 0.38192 110.29328 30.566 <2e-16 ***
## mri_age_y 0.09221 0.08992 107.78610 1.026 0.307
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## mri_age_y -0.960acg_cre_age_sex<- lmer(Cre_PCr_conc_molal ~ mri_age_y + female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_cre_age_sex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + female + (1 | subj_id)
## Data: acg_dat
##
## AIC BIC logLik deviance df.resid
## 342.9 356.5 -166.4 332.9 108
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1065 -0.3932 0.0080 0.3605 1.6020
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.6866 0.8286
## Residual 0.4871 0.6979
## Number of obs: 113, groups: subj_id, 102
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 11.65746 0.38678 111.48623 30.139 <2e-16 ***
## mri_age_y 0.08988 0.09034 107.34558 0.995 0.322
## female 0.05586 0.21445 96.67410 0.260 0.795
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y
## mri_age_y -0.927
## female -0.160 -0.100anova(acg_cre_age,acg_cre_age_sex)## Data: acg_dat
## Models:
## acg_cre_age: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_cre_age_sex: Cre_PCr_conc_molal ~ mri_age_y + female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## acg_cre_age 4 340.94 351.85 -166.47 332.94
## acg_cre_age_sex 5 342.87 356.51 -166.44 332.87 0.0678 1 0.7945acg_cre_ageXsex<- lmer(Cre_PCr_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_cre_ageXsex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + female + mri_age_y:female +
## (1 | subj_id)
## Data: acg_dat
##
## AIC BIC logLik deviance df.resid
## 343.6 360.0 -165.8 331.6 107
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.11740 -0.37592 0.04868 0.38312 1.59263
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.7009 0.8372
## Residual 0.4638 0.6810
## Number of obs: 113, groups: subj_id, 102
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.09093 0.54699 105.10134 22.104 <2e-16 ***
## mri_age_y -0.01931 0.13288 106.80701 -0.145 0.885
## female -0.77294 0.76525 111.06817 -1.010 0.315
## mri_age_y:female 0.20325 0.18016 110.02199 1.128 0.262
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y female
## mri_age_y -0.965
## female -0.715 0.689
## mri_g_y:fml 0.711 -0.738 -0.960anova(acg_cre_age,acg_cre_ageXsex)## Data: acg_dat
## Models:
## acg_cre_age: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_cre_ageXsex: Cre_PCr_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## acg_cre_age 4 340.94 351.85 -166.47 332.94
## acg_cre_ageXsex 6 343.61 359.98 -165.81 331.61 1.3273 2 0.515Cho
acg_cho_age<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_cho_age)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
## Data: acg_dat
##
## AIC BIC logLik deviance df.resid
## 78.1 89.0 -35.0 70.1 108
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1325 -0.5979 -0.1286 0.6106 2.4400
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.03204 0.1790
## Residual 0.07887 0.2808
## Number of obs: 112, groups: subj_id, 101
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.74592 0.12031 111.92574 22.82 <2e-16 ***
## mri_age_y -0.06423 0.02842 111.94035 -2.26 0.0258 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## mri_age_y -0.963acg_cho_age_sex<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_cho_age_sex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + (1 | subj_id)
## Data: acg_dat
##
## AIC BIC logLik deviance df.resid
## 78.2 91.8 -34.1 68.2 107
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.17836 -0.63945 -0.09041 0.54644 2.50032
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.03214 0.1793
## Residual 0.07694 0.2774
## Number of obs: 112, groups: subj_id, 101
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.77074 0.12063 111.75493 22.970 <2e-16 ***
## mri_age_y -0.06004 0.02834 111.86929 -2.118 0.0364 *
## female -0.09017 0.06493 103.29750 -1.389 0.1679
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y
## mri_age_y -0.931
## female -0.148 -0.107anova(acg_cho_age,acg_cho_age_sex)## Data: acg_dat
## Models:
## acg_cho_age: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_cho_age_sex: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## acg_cho_age 4 78.077 88.951 -35.039 70.077
## acg_cho_age_sex 5 78.166 91.759 -34.083 68.166 1.9111 1 0.1668acg_cho_ageXsex<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_cho_ageXsex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + mri_age_y:female +
## (1 | subj_id)
## Data: acg_dat
##
## AIC BIC logLik deviance df.resid
## 79.7 96.1 -33.9 67.7 106
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1308 -0.6109 -0.1102 0.5606 2.5638
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.03297 0.1816
## Residual 0.07579 0.2753
## Number of obs: 112, groups: subj_id, 101
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.69546 0.16726 107.68134 16.116 <2e-16 ***
## mri_age_y -0.04107 0.04072 109.53005 -1.008 0.315
## female 0.05961 0.23991 111.98464 0.248 0.804
## mri_age_y:female -0.03674 0.05661 111.93901 -0.649 0.518
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y female
## mri_age_y -0.965
## female -0.697 0.673
## mri_g_y:fml 0.694 -0.719 -0.963anova(acg_cho_age,acg_cho_ageXsex)## Data: acg_dat
## Models:
## acg_cho_age: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_cho_ageXsex: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## acg_cho_age 4 78.077 88.951 -35.039 70.077
## acg_cho_ageXsex 6 79.747 96.058 -33.874 67.747 2.3297 2 0.312#including sex in model does not significantly improve fit Ins
acg_ins_age<- lmer(Ins_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_ins_age)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
## Data: acg_dat
##
## AIC BIC logLik deviance df.resid
## 416.5 427.3 -204.2 408.5 108
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.51194 -0.47558 0.09217 0.51196 1.57868
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.9268 0.9627
## Residual 1.3749 1.1726
## Number of obs: 112, groups: subj_id, 101
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 8.4742 0.5464 111.9652 15.509 <2e-16 ***
## mri_age_y -0.2049 0.1294 111.4483 -1.583 0.116
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## mri_age_y -0.962acg_ins_age_sex<- lmer(Ins_conc_molal ~ mri_age_y + female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_ins_age_sex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + female + (1 | subj_id)
## Data: acg_dat
##
## AIC BIC logLik deviance df.resid
## 417.1 430.7 -203.5 407.1 107
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.47103 -0.48748 0.05729 0.54401 1.51171
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.8891 0.9429
## Residual 1.3819 1.1755
## Number of obs: 112, groups: subj_id, 101
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 8.5595 0.5485 111.9724 15.605 <2e-16 ***
## mri_age_y -0.1859 0.1296 111.3196 -1.435 0.154
## female -0.3489 0.2981 99.5526 -1.171 0.245
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y
## mri_age_y -0.930
## female -0.140 -0.118anova(acg_ins_age,acg_ins_age_sex)## Data: acg_dat
## Models:
## acg_ins_age: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_ins_age_sex: Ins_conc_molal ~ mri_age_y + female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## acg_ins_age 4 416.45 427.33 -204.23 408.45
## acg_ins_age_sex 5 417.09 430.69 -203.55 407.09 1.3598 1 0.2436acg_ins_ageXsex<- lmer(Ins_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_ins_ageXsex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 |
## subj_id)
## Data: acg_dat
##
## AIC BIC logLik deviance df.resid
## 418.2 434.5 -203.1 406.2 106
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.55603 -0.53421 0.07378 0.57588 1.65390
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.7615 0.8727
## Residual 1.4790 1.2162
## Number of obs: 112, groups: subj_id, 101
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 9.0615 0.7667 105.7180 11.819 <2e-16 ***
## mri_age_y -0.3138 0.1882 108.1469 -1.667 0.0983 .
## female -1.3693 1.0920 111.9981 -1.254 0.2125
## mri_age_y:female 0.2511 0.2588 111.8723 0.970 0.3340
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y female
## mri_age_y -0.965
## female -0.702 0.678
## mri_g_y:fml 0.702 -0.727 -0.963anova(acg_ins_age,acg_ins_ageXsex)## Data: acg_dat
## Models:
## acg_ins_age: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_ins_ageXsex: Ins_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## acg_ins_age 4 416.45 427.33 -204.23 408.45
## acg_ins_ageXsex 6 418.20 434.51 -203.10 406.20 2.2549 2 0.3239Glx
acg_glx_age<- lmer(Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_glx_age)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
## Data: acg_dat
##
## AIC BIC logLik deviance df.resid
## 541.1 552.0 -266.5 533.1 109
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.48535 -0.37243 -0.00651 0.43846 1.68169
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 4.706 2.169
## Residual 2.330 1.526
## Number of obs: 113, groups: subj_id, 102
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 25.7502 0.9206 106.5717 27.972 <2e-16 ***
## mri_age_y 0.2612 0.2163 103.0785 1.207 0.23
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## mri_age_y -0.959acg_glx_age_sex<- lmer(Glu_Gln_conc_molal ~ mri_age_y + female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_glx_age_sex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + female + (1 | subj_id)
## Data: acg_dat
##
## AIC BIC logLik deviance df.resid
## 543.0 556.6 -266.5 533.0 108
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.47879 -0.39442 -0.02857 0.42851 1.70104
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 4.692 2.166
## Residual 2.335 1.528
## Number of obs: 113, groups: subj_id, 102
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 25.7002 0.9334 108.5734 27.535 <2e-16 ***
## mri_age_y 0.2542 0.2173 102.5984 1.170 0.245
## female 0.1700 0.5257 98.8828 0.323 0.747
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y
## mri_age_y -0.925
## female -0.167 -0.097anova(acg_glx_age,acg_glx_age_sex)## Data: acg_dat
## Models:
## acg_glx_age: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_glx_age_sex: Glu_Gln_conc_molal ~ mri_age_y + female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## acg_glx_age 4 541.09 552.00 -266.55 533.09
## acg_glx_age_sex 5 542.99 556.62 -266.49 532.99 0.1045 1 0.7465acg_glx_ageXsex<- lmer(Glu_Gln_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_glx_ageXsex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + female + mri_age_y:female +
## (1 | subj_id)
## Data: acg_dat
##
## AIC BIC logLik deviance df.resid
## 542.5 558.9 -265.2 530.5 107
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.53333 -0.43763 0.07324 0.38166 1.63409
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 4.573 2.139
## Residual 2.297 1.515
## Number of obs: 113, groups: subj_id, 102
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 27.2150 1.3265 106.5825 20.516 <2e-16 ***
## mri_age_y -0.1271 0.3220 107.8777 -0.395 0.694
## female -2.6328 1.8384 109.4802 -1.432 0.155
## mri_age_y:female 0.6875 0.4325 108.2254 1.590 0.115
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y female
## mri_age_y -0.964
## female -0.722 0.696
## mri_g_y:fml 0.718 -0.745 -0.959anova(acg_glx_age,acg_glx_ageXsex)## Data: acg_dat
## Models:
## acg_glx_age: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_glx_ageXsex: Glu_Gln_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## acg_glx_age 4 541.09 552.00 -266.55 533.09
## acg_glx_ageXsex 6 542.49 558.85 -265.25 530.49 2.603 2 0.2721Plot results ACG
Create scatter plot with regression line based on model predictions NAA
#pull intercept and slope from model for regression line
fixef.acg_naa<-fixef(acg_naa_age)
#plot
plot_acg_naa_age<-ggplot(acg_dat, aes(x=mri_age_y, y=NAA_conc_molal,color=as.factor(female)))
plot_acg_naa_age +
geom_point() +
geom_line(aes(group = as.factor(subj_id))) +
labs(title = "Anterior Cingulate", x = "Age (years)", y = "tNAA (molal)") +
scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
geom_abline(aes(intercept = fixef.acg_naa[[1]], slope= fixef.acg_naa[[2]])) +
theme_classic()
ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/acg_naa_age.png")## Saving 7 x 5 in imageCho
#pull intercept and slope from model for regression line
fixef.acg_cho<-fixef(acg_cho_age)
#plot
plot_acg_cho_age<-ggplot(acg_dat, aes(x=mri_age_y, y=Cho_GPC_PCh_conc_molal,color=as.factor(female)))
plot_acg_cho_age +
geom_point() +
geom_line(aes(group = as.factor(subj_id))) +
labs(title = "Anterior Cingulate", x = "Age (years)", y = "tCho (molal)") +
scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
geom_abline(aes(intercept = fixef.acg_cho[[1]], slope= fixef.acg_cho[[2]])) +
theme_classic()## Warning: Removed 1 rows containing missing values (geom_point).## Warning: Removed 1 row(s) containing missing values (geom_path).
ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/acg_cho_age.png")## Saving 7 x 5 in image## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).Ins
#pull intercept and slope from model for regression line
fixef.acg_ins<-fixef(acg_ins_age)
#plot
plot_acg_ins_age<-ggplot(acg_dat, aes(x=mri_age_y, y=Ins_conc_molal,color=as.factor(female)))
plot_acg_ins_age +
geom_point() +
geom_line(aes(group = as.factor(subj_id))) +
labs(title = "Anterior Cingulate", x = "Age (years)", y = "Ins (molal)") +
scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
geom_abline(aes(intercept = fixef.acg_ins[[1]], slope= fixef.acg_ins[[2]])) +
theme_classic()## Warning: Removed 1 rows containing missing values (geom_point).## Warning: Removed 1 row(s) containing missing values (geom_path).
ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/acg_ins_age.png")## Saving 7 x 5 in image## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).LME modelling LAG
run mixed-effects models for lag with participant as random effect, age, sex and age x sex interaction as fixed effects NAA, substantial longitudinal data so we include random slope run on lag_data_molal_outliers_removed_n320.csv
#basic model effect of age on intercept
lag_naa_age<- lmer(NAA_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_naa_age)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## Data: lag_dat
##
## AIC BIC logLik deviance df.resid
## 896 911 -444 888 314
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4313 -0.5792 -0.0185 0.5440 3.6263
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.1641 0.4051
## Residual 0.8272 0.9095
## Number of obs: 318, groups: subj_id, 95
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 13.82375 0.18588 309.72511 74.369 < 2e-16 ***
## mri_age_y 0.09699 0.02902 316.51634 3.343 0.00093 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## mri_age_y -0.929lag_naa_age_sex<- lmer(NAA_conc_molal ~ mri_age_y + female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_naa_age_sex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + female + (1 | subj_id)
## Data: lag_dat
##
## AIC BIC logLik deviance df.resid
## 898.0 916.8 -444.0 888.0 313
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4286 -0.5791 -0.0185 0.5452 3.6231
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.1637 0.4046
## Residual 0.8275 0.9097
## Number of obs: 318, groups: subj_id, 95
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 13.827819 0.194854 290.918032 70.965 < 2e-16 ***
## mri_age_y 0.097086 0.029038 316.589446 3.343 0.000927 ***
## female -0.009878 0.137557 93.706753 -0.072 0.942903
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y
## mri_age_y -0.874
## female -0.300 -0.039anova(lag_naa_age,lag_naa_age_sex)## Data: lag_dat
## Models:
## lag_naa_age: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_naa_age_sex: NAA_conc_molal ~ mri_age_y + female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## lag_naa_age 4 895.96 911.01 -443.98 887.96
## lag_naa_age_sex 5 897.95 916.76 -443.98 887.95 0.0051 1 0.943lag_naa_ageXsex<- lmer(NAA_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_naa_ageXsex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 |
## subj_id)
## Data: lag_dat
##
## AIC BIC logLik deviance df.resid
## 899.9 922.5 -444.0 887.9 312
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4372 -0.5859 -0.0225 0.5510 3.6185
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.1645 0.4056
## Residual 0.8268 0.9093
## Number of obs: 318, groups: subj_id, 95
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 13.79473 0.24812 314.11810 55.598 < 2e-16 ***
## mri_age_y 0.10272 0.03910 313.83715 2.627 0.00903 **
## female 0.06577 0.37459 308.18552 0.176 0.86075
## mri_age_y:female -0.01265 0.05838 316.98288 -0.217 0.82864
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y female
## mri_age_y -0.924
## female -0.662 0.612
## mri_g_y:fml 0.619 -0.670 -0.930anova(lag_naa_age,lag_naa_ageXsex)## Data: lag_dat
## Models:
## lag_naa_age: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_naa_ageXsex: NAA_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## lag_naa_age 4 895.96 911.01 -443.98 887.96
## lag_naa_ageXsex 6 899.91 922.48 -443.95 887.91 0.0518 2 0.9744#including sex in model does not significantly improve fit
#test quadratic effect of age
#create age^2 variable
lag_dat$mri_age_y2 <- (lag_dat$mri_age_y)^2
lag_naa_age2<- lmer(NAA_conc_molal ~ mri_age_y + mri_age_y2 + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_naa_age2) #age^2 p=.0447## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + mri_age_y2 + (1 | subj_id)
## Data: lag_dat
##
## AIC BIC logLik deviance df.resid
## 893.9 912.7 -442.0 883.9 313
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5451 -0.5819 -0.0270 0.4941 3.5330
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.1700 0.4123
## Residual 0.8119 0.9011
## Number of obs: 318, groups: subj_id, 95
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.70181 0.58588 310.95949 21.680 <2e-16 ***
## mri_age_y 0.47130 0.18749 303.10334 2.514 0.0125 *
## mri_age_y2 -0.02836 0.01402 297.51463 -2.023 0.0440 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y
## mri_age_y -0.983
## mri_age_y2 0.949 -0.988anova(lag_naa_age, lag_naa_age2) #fit improved p=0.44## Data: lag_dat
## Models:
## lag_naa_age: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_naa_age2: NAA_conc_molal ~ mri_age_y + mri_age_y2 + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## lag_naa_age 4 895.96 911.01 -443.98 887.96
## lag_naa_age2 5 893.91 912.72 -441.96 883.91 4.0455 1 0.04429 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Cre
lag_cre_age<- lmer(Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cre_age)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
## Data: lag_dat
##
## AIC BIC logLik deviance df.resid
## 848.2 863.3 -420.1 840.2 316
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9764 -0.5307 -0.1261 0.5525 4.1548
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.08477 0.2912
## Residual 0.73683 0.8584
## Number of obs: 320, groups: subj_id, 95
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 10.49912 0.16805 308.49608 62.476 <2e-16 ***
## mri_age_y -0.04183 0.02646 319.87112 -1.581 0.115
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## mri_age_y -0.937lag_cre_age_sex<- lmer(Cre_PCr_conc_molal ~ mri_age_y + female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cre_age_sex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + female + (1 | subj_id)
## Data: lag_dat
##
## AIC BIC logLik deviance df.resid
## 849.7 868.6 -419.9 839.7 315
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9528 -0.5658 -0.1208 0.5434 4.1925
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.08478 0.2912
## Residual 0.73557 0.8577
## Number of obs: 320, groups: subj_id, 95
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 10.46482 0.17478 289.38001 59.873 <2e-16 ***
## mri_age_y -0.04267 0.02646 319.88902 -1.612 0.108
## female 0.08318 0.11739 84.73092 0.709 0.481
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y
## mri_age_y -0.887
## female -0.277 -0.044anova(lag_cre_age,lag_cre_age_sex)## Data: lag_dat
## Models:
## lag_cre_age: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_cre_age_sex: Cre_PCr_conc_molal ~ mri_age_y + female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## lag_cre_age 4 848.23 863.31 -420.12 840.23
## lag_cre_age_sex 5 849.73 868.57 -419.87 839.73 0.5017 1 0.4788lag_cre_ageXsex<- lmer(Cre_PCr_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cre_ageXsex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + female + mri_age_y:female +
## (1 | subj_id)
## Data: lag_dat
##
## AIC BIC logLik deviance df.resid
## 850.0 872.6 -419.0 838.0 314
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9579 -0.5653 -0.0910 0.5718 4.2323
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.08177 0.2860
## Residual 0.73339 0.8564
## Number of obs: 320, groups: subj_id, 95
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 10.65304 0.22484 314.11103 47.381 <2e-16 ***
## mri_age_y -0.07476 0.03585 318.90184 -2.086 0.0378 *
## female -0.33761 0.33688 307.33953 -1.002 0.3171
## mri_age_y:female 0.07047 0.05297 319.97938 1.330 0.1844
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y female
## mri_age_y -0.935
## female -0.667 0.624
## mri_g_y:fml 0.632 -0.677 -0.938anova(lag_cre_age,lag_cre_ageXsex)## Data: lag_dat
## Models:
## lag_cre_age: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_cre_ageXsex: Cre_PCr_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## lag_cre_age 4 848.23 863.31 -420.12 840.23
## lag_cre_ageXsex 6 849.97 872.58 -418.98 837.97 2.263 2 0.3226Cho
lag_cho_age<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cho_age)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
## Data: lag_dat
##
## AIC BIC logLik deviance df.resid
## 51.3 66.4 -21.7 43.3 316
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6539 -0.5810 -0.0562 0.4659 4.1202
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.0149 0.1221
## Residual 0.0561 0.2368
## Number of obs: 320, groups: subj_id, 95
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.346670 0.049038 310.144986 47.854 < 2e-16 ***
## mri_age_y -0.023315 0.007581 315.773363 -3.075 0.00229 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## mri_age_y -0.921lag_cho_age_sex<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cho_age_sex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + (1 | subj_id)
## Data: lag_dat
##
## AIC BIC logLik deviance df.resid
## 51.6 70.5 -20.8 41.6 315
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6317 -0.5893 -0.0587 0.4688 4.1717
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.01469 0.1212
## Residual 0.05586 0.2364
## Number of obs: 320, groups: subj_id, 95
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.36786 0.05147 283.09492 46.006 < 2e-16 ***
## mri_age_y -0.02286 0.00757 315.96555 -3.019 0.00274 **
## female -0.05024 0.03808 80.73098 -1.319 0.19083
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y
## mri_age_y -0.860
## female -0.312 -0.046anova(lag_cho_age,lag_cho_age_sex)## Data: lag_dat
## Models:
## lag_cho_age: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_cho_age_sex: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## lag_cho_age 4 51.345 66.419 -21.673 43.345
## lag_cho_age_sex 5 51.610 70.452 -20.805 41.610 1.735 1 0.1878lag_cho_ageXsex<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cho_ageXsex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + mri_age_y:female +
## (1 | subj_id)
## Data: lag_dat
##
## AIC BIC logLik deviance df.resid
## 49.9 72.5 -18.9 37.9 314
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6639 -0.5411 -0.0413 0.4642 4.0526
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.01452 0.1205
## Residual 0.05521 0.2350
## Number of obs: 320, groups: subj_id, 95
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.28998 0.06505 314.68723 35.206 <2e-16 ***
## mri_age_y -0.00954 0.01019 312.12139 -0.936 0.3497
## female 0.12488 0.09791 309.31059 1.275 0.2031
## mri_age_y:female -0.02931 0.01511 316.65469 -1.939 0.0533 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y female
## mri_age_y -0.916
## female -0.664 0.609
## mri_g_y:fml 0.617 -0.674 -0.922anova(lag_cho_age,lag_cho_ageXsex)## Data: lag_dat
## Models:
## lag_cho_age: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_cho_ageXsex: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## lag_cho_age 4 51.345 66.419 -21.673 43.345
## lag_cho_ageXsex 6 49.871 72.481 -18.936 37.871 5.4741 2 0.06476 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#including sex in model does not significantly improve fit Ins
lag_ins_age<- lmer(Ins_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_ins_age)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
## Data: lag_dat
##
## AIC BIC logLik deviance df.resid
## 1013.3 1028.3 -502.6 1005.3 315
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3984 -0.5594 -0.0111 0.5911 4.6516
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.1768 0.4205
## Residual 1.2227 1.1057
## Number of obs: 319, groups: subj_id, 95
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 6.28473 0.21938 307.28604 28.648 <2e-16 ***
## mri_age_y 0.01632 0.03444 318.45543 0.474 0.636
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## mri_age_y -0.934lag_ins_age_sex<- lmer(Ins_conc_molal ~ mri_age_y + female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_ins_age_sex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + female + (1 | subj_id)
## Data: lag_dat
##
## AIC BIC logLik deviance df.resid
## 1015.3 1034.1 -502.6 1005.3 314
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4003 -0.5588 -0.0120 0.5901 4.6506
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.1765 0.4201
## Residual 1.2229 1.1058
## Number of obs: 319, groups: subj_id, 95
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 6.286355 0.228788 284.644902 27.477 <2e-16 ***
## mri_age_y 0.016356 0.034473 318.479871 0.474 0.635
## female -0.003854 0.157311 78.609245 -0.024 0.981
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y
## mri_age_y -0.881
## female -0.284 -0.046anova(lag_ins_age,lag_ins_age_sex)## Data: lag_dat
## Models:
## lag_ins_age: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_ins_age_sex: Ins_conc_molal ~ mri_age_y + female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## lag_ins_age 4 1013.3 1028.3 -502.64 1005.3
## lag_ins_age_sex 5 1015.3 1034.1 -502.64 1005.3 6e-04 1 0.9807lag_ins_ageXsex<- lmer(Ins_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_ins_ageXsex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 |
## subj_id)
## Data: lag_dat
##
## AIC BIC logLik deviance df.resid
## 1017.1 1039.7 -502.6 1005.1 313
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4248 -0.5418 -0.0121 0.5899 4.6348
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.1767 0.4204
## Residual 1.2221 1.1055
## Number of obs: 319, groups: subj_id, 95
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 6.361751 0.294417 312.731852 21.608 <2e-16 ***
## mri_age_y 0.003458 0.046822 316.620308 0.074 0.941
## female -0.171702 0.441468 305.644875 -0.389 0.698
## mri_age_y:female 0.028161 0.069170 318.668940 0.407 0.684
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y female
## mri_age_y -0.930
## female -0.667 0.620
## mri_g_y:fml 0.630 -0.677 -0.934anova(lag_ins_age,lag_ins_ageXsex)## Data: lag_dat
## Models:
## lag_ins_age: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_ins_ageXsex: Ins_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## lag_ins_age 4 1013.3 1028.3 -502.64 1005.3
## lag_ins_ageXsex 6 1017.1 1039.7 -502.56 1005.1 0.1663 2 0.9202Glx
lag_glx_age<- lmer(Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)## boundary (singular) fit: see ?isSingularisSingular(lag_glx_age) # TRUE: model is almost/near singular (parameters are on the boundary of the feasible parameter space: random effects variance = 0)## [1] TRUEsummary(lag_glx_age)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
## Data: lag_dat
##
## AIC BIC logLik deviance df.resid
## 1396.1 1411.1 -694.0 1388.1 315
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6700 -0.6848 0.0127 0.6729 3.5338
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.000 0.000
## Residual 4.542 2.131
## Number of obs: 319, groups: subj_id, 95
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 23.23676 0.39326 319.00000 59.088 < 2e-16 ***
## mri_age_y -0.25961 0.06309 319.00000 -4.115 4.93e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## mri_age_y -0.953
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularlag_glx_age_sex<- lmer(Glu_Gln_conc_molal ~ mri_age_y + female + (1 | subj_id), data=lag_dat, REML=FALSE)## boundary (singular) fit: see ?isSingularsummary(lag_glx_age_sex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + female + (1 | subj_id)
## Data: lag_dat
##
## AIC BIC logLik deviance df.resid
## 1396.7 1415.5 -693.4 1386.7 314
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7491 -0.6907 0.0192 0.7055 3.6028
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.000 0.000
## Residual 4.523 2.127
## Number of obs: 319, groups: subj_id, 95
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 23.12126 0.40478 319.00000 57.120 < 2e-16 ***
## mri_age_y -0.26202 0.06299 319.00000 -4.160 4.1e-05 ***
## female 0.27789 0.23880 319.00000 1.164 0.245
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y
## mri_age_y -0.915
## female -0.245 -0.033
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularanova(lag_glx_age,lag_glx_age_sex)## Data: lag_dat
## Models:
## lag_glx_age: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_glx_age_sex: Glu_Gln_conc_molal ~ mri_age_y + female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## lag_glx_age 4 1396.1 1411.1 -694.04 1388.1
## lag_glx_age_sex 5 1396.7 1415.5 -693.36 1386.7 1.3513 1 0.245lag_glx_ageXsex<- lmer(Glu_Gln_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=lag_dat, REML=FALSE)## boundary (singular) fit: see ?isSingularsummary(lag_glx_ageXsex)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + female + mri_age_y:female +
## (1 | subj_id)
## Data: lag_dat
##
## AIC BIC logLik deviance df.resid
## 1397.8 1420.4 -692.9 1385.8 313
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7059 -0.6942 -0.0013 0.6785 3.6197
##
## Random effects:
## Groups Name Variance Std.Dev.
## subj_id (Intercept) 0.00 0.000
## Residual 4.51 2.124
## Number of obs: 319, groups: subj_id, 95
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 23.44862 0.52588 319.00000 44.589 < 2e-16 ***
## mri_age_y -0.31768 0.08502 319.00000 -3.737 0.000221 ***
## female -0.45362 0.78868 319.00000 -0.575 0.565583
## mri_age_y:female 0.12295 0.12636 319.00000 0.973 0.331265
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mr_g_y female
## mri_age_y -0.951
## female -0.667 0.634
## mri_g_y:fml 0.640 -0.673 -0.953
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingularanova(lag_glx_age,lag_glx_ageXsex)## Data: lag_dat
## Models:
## lag_glx_age: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_glx_ageXsex: Glu_Gln_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## lag_glx_age 4 1396.1 1411.1 -694.04 1388.1
## lag_glx_ageXsex 6 1397.8 1420.4 -692.89 1385.8 2.2967 2 0.3172Plot results LAG
Create scatter plot with regression line based on model predictions NAA
#pull intercept and slope from model for regression line
fixef.lag_naa<-fixef(lag_naa_age)
#plot
plot_lag_naa_age<-ggplot(lag_dat, aes(x=mri_age_y, y=NAA_conc_molal,color=as.factor(female)))
plot_lag_naa_age +
geom_point() +
geom_line(aes(group = as.factor(subj_id))) +
labs(title = "Left Temporo-Parietal", x = "Age (years)", y = "tNAA (molal)") +
scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
geom_abline(aes(intercept = fixef.lag_naa[[1]], slope= fixef.lag_naa[[2]])) +
theme_classic()## Warning: Removed 3 rows containing missing values (geom_point).## Warning: Removed 3 row(s) containing missing values (geom_path).
ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/lag_naa_age.png")## Saving 7 x 5 in image## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 3 row(s) containing missing values (geom_path).Cho
#pull intercept and slope from model for regression line
fixef.lag_cho<-fixef(lag_cho_age)
#plot
plot_lag_cho_age<-ggplot(lag_dat, aes(x=mri_age_y, y=Cho_GPC_PCh_conc_molal,color=as.factor(female)))
plot_lag_cho_age +
geom_point() +
geom_line(aes(group = as.factor(subj_id))) +
labs(title = "Left Temporo-Parietal", x = "Age (years)", y = "tCho (molal)") +
scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
geom_abline(aes(intercept = fixef.lag_cho[[1]], slope= fixef.lag_cho[[2]])) +
theme_classic()## Warning: Removed 1 rows containing missing values (geom_point).## Warning: Removed 1 row(s) containing missing values (geom_path).
ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/lag_cho_age.png")## Saving 7 x 5 in image## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).Glx
#pull intercept and slope from model for regression line
fixef.lag_glx<-fixef(lag_glx_age)
#plot
plot_lag_glx_age<-ggplot(lag_dat, aes(x=mri_age_y, y=Glu_Gln_conc_molal,color=as.factor(female)))
plot_lag_glx_age +
geom_point() +
geom_line(aes(group = as.factor(subj_id))) +
labs(title = "Left Temporo-Parietal", x = "Age (years)", y = "Glx (molal)") +
scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
geom_abline(aes(intercept = fixef.lag_glx[[1]], slope= fixef.lag_glx[[2]])) +
theme_classic()## Warning: Removed 2 rows containing missing values (geom_point).## Warning: Removed 2 row(s) containing missing values (geom_path).
ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/lag_glx_age.png")## Saving 7 x 5 in image## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 2 row(s) containing missing values (geom_path).##Correct for multiple comparisons Use FDR method (Benjamini & Hochberg) to correct for mutliple comparisons (5 metabolites x 2 voxels = 10 models)
#list models (voxel/metabolite pair)
models<- c("acg_naa", "acg_cr", "acg_cho", "acg_ins", "acg_glx", "lag_naa", "lag_cr", "lag_cho", "lag_ins", "lag_glx")
#list p-values for each model, matching order of models
p<-c(.0217, .307, .0258, .0494, .23, .00093, .115, .00229, .636, 0)
fdr_table<-data.frame(models, p)
fdr_table$p_fdr<-p.adjust(p, method = "fdr", n= length(p))
fdr_table## models p p_fdr
## 1 acg_naa 0.02170 0.051600000
## 2 acg_cr 0.30700 0.341111111
## 3 acg_cho 0.02580 0.051600000
## 4 acg_ins 0.04940 0.082333333
## 5 acg_glx 0.23000 0.287500000
## 6 lag_naa 0.00093 0.004650000
## 7 lag_cr 0.11500 0.164285714
## 8 lag_cho 0.00229 0.007633333
## 9 lag_ins 0.63600 0.636000000
## 10 lag_glx 0.00000 0.000000000#run FDR without LAG-Glx in case it's invalid:
p2<-p<-c(.0217, .307, .0258, .0494, .23, .00093, .115, .00229, .636)
p.adjust(p2, method = "fdr", n= length(p2))## [1] 0.0580500 0.3453750 0.0580500 0.0889200 0.2957143 0.0083700 0.1725000
## [8] 0.0103050 0.6360000##Combine plots into 1 figure Simplify plots to legend & region don’t repeat in each panel
library(patchwork)
p1_acg_naa_age<-plot_acg_naa_age +
geom_point() +
geom_line(aes(group = as.factor(subj_id))) +
labs(title = "Anterior Cingulate", x = "Age (years)", y = "tNAA (molal)") +
scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
geom_abline(aes(intercept = fixef.acg_naa[[1]], slope= fixef.acg_naa[[2]])) +
theme_classic() +
theme(plot.title= element_text(hjust = 0.5))
p4_lag_naa_age<-plot_lag_naa_age +
geom_point() +
geom_line(aes(group = as.factor(subj_id))) +
labs(title = "Left Temporo-Parietal", x = "Age (years)", y = "tNAA (molal)") +
scale_color_discrete(guide="legend", name = "Sex", labels = c("Male", "Female")) +
theme_classic() +
#theme(legend.position = c(0.7, 0.95), legend.background = element_rect(fill='transparent')) +
geom_abline(aes(intercept = fixef.lag_naa[[1]], slope= fixef.lag_naa[[2]])) +
theme(plot.title= element_text(hjust = 0.5))
p2_acg_cho_age<-plot_acg_cho_age +
geom_point() +
geom_line(aes(group = as.factor(subj_id))) +
labs(x = "Age (years)", y = "tCho (molal)") +
scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
geom_abline(aes(intercept = fixef.acg_cho[[1]], slope= fixef.acg_cho[[2]])) +
theme_classic()
p5_lag_cho_age<-plot_lag_cho_age +
geom_point() +
geom_line(aes(group = as.factor(subj_id))) +
labs(x = "Age (years)", y = "tCho (molal)") +
scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
geom_abline(aes(intercept = fixef.lag_cho[[1]], slope= fixef.lag_cho[[2]])) +
theme_classic()
p3_acg_ins_age<-plot_acg_ins_age +
geom_point() +
geom_line(aes(group = as.factor(subj_id))) +
labs(x = "Age (years)", y = "Ins (molal)") +
scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
geom_abline(aes(intercept = fixef.acg_ins[[1]], slope= fixef.acg_ins[[2]])) +
theme_classic()
p6_lag_glx_age<-plot_lag_glx_age +
geom_point() +
geom_line(aes(group = as.factor(subj_id))) +
labs(x = "Age (years)", y = "Glx (molal)") +
scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
geom_abline(aes(intercept = fixef.lag_glx[[1]], slope= fixef.lag_glx[[2]])) +
theme_classic()
results_fig<-p1_acg_naa_age + p4_lag_naa_age + p2_acg_cho_age + p5_lag_cho_age + p3_acg_ins_age + p6_lag_glx_age + plot_layout(nrow=3, guides ='collect')
ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/results_fig_grid.png", results_fig, width = 2000, units = "px")## Saving 2000 x 1500 px image## Warning: Removed 3 rows containing missing values (geom_point).## Warning: Removed 3 row(s) containing missing values (geom_path).## Warning: Removed 1 rows containing missing values (geom_point).## Warning: Removed 1 row(s) containing missing values (geom_path).## Warning: Removed 1 rows containing missing values (geom_point).## Warning: Removed 1 row(s) containing missing values (geom_path).## Warning: Removed 1 rows containing missing values (geom_point).## Warning: Removed 1 row(s) containing missing values (geom_path).## Warning: Removed 2 rows containing missing values (geom_point).## Warning: Removed 2 row(s) containing missing values (geom_path).